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Abstract 

A  research  approach  is  presented  that  extends  the  separate  solntion  methods  of 
stochastic  and  mnlti-objective  optimization  problems  to  one  that  wonld  solve  prob¬ 
lems  having  both  characteristics.  Snch  problems  are  typically  enconntered  when 
one  desires  to  optimize  systems  with  mnltiple,  often  competing,  objectives  that  do 
not  have  a  closed  form  representation  and  mnst  be  estimated,  e.g.,  via  simnlation. 
First,  the  class  of  mesh  adaptive  direct  search  (MADS)  algorithms  for  nonlinearly 
constrained  optimization  is  extended  to  mixed  variable  problems,  and  convergence 
to  appropriately  dehned  stationary  points  is  proved.  The  resnlting  algorithm,  MV- 
MADS,  is  then  extended  to  stochastic  problems  (MVMADS-RS),  via  a  ranking  and 
selection  procednre.  Finally,  a  two-stage  method  is  developed  that  combines  the 
generalized  pattern  search/ranking  and  selection  (MGPS-RS)  algorithms  for  single¬ 
objective,  mixed  variable,  stochastic  problems  with  a  mnlti-objective  approach  that 
makes  nse  of  interactive  techniqnes  for  the  specihcation  of  aspiration  and  reservation 
levels,  scalarization  fnnctions,  and  mnlti-objective  ranking  and  selection.  This  com¬ 
bination  is  devised  specihcally  so  as  to  keep  the  desirable  convergence  properties  of 
MGPS-RS  and  MVMADS-RS,  while  extending  to  the  mnlti-objective  case.  A  conver¬ 
gence  analysis  for  the  general  class  of  algorithms  establishes  almost  snre  convergence 
of  an  iteration  snbseqnence  to  stationary  points  appropriately  dehned  in  the  mixed- 
variable  domain.  Seven  specihc  instances  of  the  new  algorithm  are  implemented 
and  tested  on  11  mnlti-objective  test  problems  from  the  literatnre  and  an  engineering 
design  problem. 
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Search  Techniques  eor  Multi- Objective  Optimization  oe 
Mixed- Variable  Systems  Having  Stochastic  Responses 

I.  Introduction 

1.1  Problem  Setting 

With  the  advent  of  advanced  nnmerical  techniqnes  for  analyzing  complex  engi¬ 
neering  problems,  engineers  are  seeking  to  integrate  snch  methods  into  smart  design 
processes.  As  a  resnlt,  optimization  techniqnes  have  become  increasingly  important 
in  engineering  design  problems,  e.g.,  in  the  area  of  conceptnal  aircraft  design  [99]. 
However,  complications  arise  when  applying  traditional  optimization  techniqnes  to 
engineering  design  problems  like  aircraft  design.  Snch  design  problems  typically 
contain  mnltiple,  often  competing,  objectives  [35,77,82,90].  Additionally,  these 
objectives  are  often  snbject  to  measnrement  error  or  mnst  be  estimated  by  complex 
simnlations  [35,36,81,91].  Finally,  engineering  design  problems  often  contain  discrete 
or  categorical  design  variables  [82].  Thns,  mnlti-objective  and  stochastic  optimiza¬ 
tion  techniqnes,  applicable  to  both  continnons  and  discrete  decision  variables,  mnst 
be  considered.  However,  signihcant  challenges  exist  for  these  types  of  optimization 
problems. 

1.1.1  Modeling  Uncertainty.  Given  a  general  constrained  nonlinear  opti¬ 
mization  problem,  formnlated  as 


minF(a;) 

snbject  to 

gi{x)<0,  i  =  l,2,...,M, 

where  x  G  represents  the  controllable  design  variables,  F  :  R”  —S’  R,  and  the 
constraints  gi  :  R”  — ?•  R,  ?  =  1,  2, . . . ,  M,  dehne  the  feasible  region  Q  C  R",  a  myriad 
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of  archetypal  solution  methods  are  available  to  the  analyst  {e.g.,  Newton’s  method, 
steepest  descent,  etc.)  [162],  But  if  one  or  more  elements  of  the  optimization  problem 
contain  a  random  component,  what  changes  to  classical  optimization  are  required  in 
order  to  allow  reasonably  efficient  optimization  of  this  random  system?  The  answer 
to  this  question  involves  the  class  of  solution  methods  called  stochastic  optimization 
which,  as  the  name  implies,  deals  with  systems  in  which  one  or  more  of  the  problem 
components — objective  function(s)  and/or  constraints — contain  a  random  element. 
Stochastic  optimization  is  related  to  traditional  statistical  modeling,  in  that  one  seeks 
to  model  a  function  by  repeated  sampling  of  a  population;  however,  in  this  case,  the 
function  is  dependent  not  only  on  a  random  element,  but  also  on  controllable  design 
variables. 

A  general  stochastic  optimization  problem  can  be  written  as 

minF(a;,cu)  (1.2a) 

subject  to 

gi{x,uj)<Q,  i  =  l,2,...,M  (1.2b) 

where  x  G  and  uj  G  represent  the  controllable  design  variables  and  random 
(uncontrollable  environmental)  variables,  respectively,  and  the  constraint  set  gi{x,u!), 
i  =  1,  2, . . . ,  M  dehnes  some  feasible  region  hi  C  (R"-i  x  Thus,  the  task  becomes 

hnding  x  so  that  Equation  (1.2a)  is  minimized  in  a  certain  sense  over  all  feasible  x  and 
all  possible  values  of  cu.  Notions  of  feasibility  and  optimality  for  stochastic  systems 
are  highly  dependent  on  the  specihc  problem  under  study  and  must  be  precisely 
dehned  [48]. 

This  research  focuses  on  problems  in  which  all  constraints  are  deterministic  and 
the  objective  function  cannot  be  explicitly  evaluated  and  must  be  estimated  through 
some  kind  of  simulation.  Here  simulation  refers  to  a  generic  numerical  method  by 
which  input  (control)  variables  are  used  to  produce  an  output  measure  of  interest 
(response)  [53,138].  Therefore,  in  this  simulation-based  optimization,  the  observed 
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system  response  is  a  function  of  both  the  design  variables  and  the  random  error 
associated  with  the  simulation  model. 

For  simulation-based  optimization,  the  general  stochastic  objective  function 
given  in  (1.2a)  is  typically  replaced  by  its  mathematical  expectation  E[-]  [48].  With 
this  convention  and  the  assumption  that  the  observed  response  is  an  unbiased  approx¬ 
imation  of  the  true  system  response,  the  observed  response  F  can  be  represented  by 
F(x,uj)  =  f(x)  +  ei_j{x)  where  /  is  the  deterministic,  “true”  objective  function,  and 
:  R"'!  — >  R  is  the  random  error  function  associated  with  the  simulation  (he.,  the 
random  element  of  F  due  to  the  uncontrollable  environmental  variables  cu)  such  that 
F[et^{x)\  =  0.  Thus,  the  problem  is  reformulated  as 

minF(a;,a;)  ~  F[F(x)]  =  F[f{x)  +£a;(a:)] 

subject  to 

gi{x)  <0,  i  =  1,2, . . .  ,M. 

1.1.2  Optimizing  Multiple  Objectives.  In  addition  to  the  complexity  of 
design  optimization  problems  due  to  the  stochastic  element,  often  there  exists  no 
single  criterion  for  choosing  the  best  solution.  In  fact,  even  the  notion  of  “best”  can 
be  unclear  when  multiple  objectives  are  present.  In  many  cases,  it  can  be  shown  that 
improvement  to  one  objective  actually  degrades  the  performance  of  another.  Such 
objectives  are  called  competing  objectives  and  are  the  motivation  for  the  study  of 
multi-objective  optimization  as  a  separate  discipline;  if  the  objectives  did  not  compete, 
they  could  be  collapsed  into  a  single  objective  and  classic  solution  techniques  would 
apply  [I].  The  multi-objective  optimization  problem, 

mini?[F(a;)]  =  F[f{x)  +£^^(0;)]  (I.4a) 

subject  to 

a;  G  =  {a;  e  R"^  :  gi{x)  <0,  i  =  1,2, . . . ,  M],  (I.4b) 
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where  F  :  R”  — ?•  is  that  of  hnding  a  solution  a;*  G  hi  that  optimizes  the  set  of 

objectives  F  =  (Fi,  F2, . . . ,  Fj)  in  the  sense  that  no  other  point  y  E  Vi  yields  a  better 
function  value  in  all  the  objectives  [83].  The  point  x*  is  said  to  be  non- dominated, 
efficient  or  optimal  in  the  Pareto  sense  [54],  Further  dehnitions  associated  with 
Pareto  optimality  follow. 

1.1. 2.1  Definitions  Associated  with  Pareto  Optimality.  When  optimiz¬ 
ing  over  multiple  objectives,  the  space  containing  the  decision  variables  is  referred  to 
as  the  decision  space,  and  the  space  containing  their  associated  images  (i.e.,  objective 
function  values)  is  referred  to  as  the  objective  space.  Two  particular  points  in  the 
objective  space,  which  are  often  referenced  in  multi-objective  optimization  methods, 
are  the  utopia  point  (or  ideal  point)  and  nadir  point.  These  are  formally  dehned  in 
Definitions  1.1.1  and  1.1.2,  respectively  [46],  and  are  depicted  in  Figure  1.1  by  points 
U  and  N,  respectively. 

Definition  1.1.1.  The  point  y^  =  [y^ , . . .  ,yj)  given  by  y’fi  :=  minFfc(a;)  is  called  the 
utopia  or  ideal  point  of  the  multi-objective  optimization  problem  min(Fi(a;), . . . ,  Fj{x)). 

Definition  1.1.2.  The  point  y^  =  [y^ , . . .  ,yj)  given  by  y^  :=  maxFfc(a;)  is  called 

x£Cl 

the  nadir  point  of  the  multi-objective  optimization  problem  min(Fi(a;), . . . ,  Fj{x)). 

x£fl 


Figure  1.1: 


Graphical  Depiction  of  the  Objective  Space  (Convex  Pareto  Front) 
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There  are  several  equivalent  definitions  of  Pareto  optimal  solutions — also  called  non- 
dominated  or  efficient  solutions — in  the  literature.  For  the  purpose  of  this  research, 
the  definition  presented  as  Definition  1.1.3  is  used  [46]. 

Definition  1.1.3.  A  solution  x*  E  ft  to  the  multi-objective  optimization  problem, 
given  in  Equations  (Lja-l-jb),  is  said  to  be  Pareto  optimal  if  there  is  no  x  E  ft 
such  that  Fk{x)  <  Fk{x*)  for  all  k  =  1,  2, . . . ,  J  with  Fi{x)  <  Fi{x*)  for  some  i  E 
{1,2,...,J}. 

Definition  1.1.4.  A  solution  x*  E  ft  to  the  multi-objective  optimization  problem, 
given  in  Equations  (l.ja-Fjb),  is  said  to  be  a  Pareto  stationary  point  or  to  satisfy 
first-order  necessary  conditions  for  Pareto  optimality  if  there  exists  no  feasible  direc¬ 
tion  d  E  such  that  VFk{x*ffid  <  0  for  all  k  =  1,  2, . . . ,  J  and  VFi{x*ffid  <  0  for 
some  i  E  {1,2, . . . ,  J} . 

In  this  research,  the  set  of  Pareto  optimal  solutions  is  referred  to  as  the  Pareto 
optimal  set  or  simply  the  Pareto  set.  The  image  of  the  Pareto  set  is  referred  to  as 
the  Pareto  Frontier  or  Pareto  Front.  A  graphical  example  of  a  convex  Pareto  front  is 
shown  in  Figure  1.1;  a  non-convex  Pareto  front  is  shown  in  Figure  1.2.  If  the  Pareto 
set  (or  corresponding  Pareto  front)  results  from  a  solution  algorithm  and  is  not  exact, 
it  is  referred  to  as  the  approximate  (or  experimental)  Pareto  set  or  approximate  (or 
experimental)  Pareto  front,  respectively. 

1.1.3  Optimizing  Over  a  Discrete  Decision  Space.  The  complexity  of  these 
problems  is  further  increased  when  the  decision  space  includes  variables  that  are  either 
discrete  {e.g.,  integer- valued)  or  categorical.  Categorical  variables  are  those  which 
can  only  take  on  values  from  a  predetermined  list,  and  may  not  even  have  an  ordinal 
relationship  to  each  other.  However,  since  categorical  variables  can  be  mapped  easily 
to  discrete-numeric  values,  these  two  types  of  variables  will  be  grouped  together  and 
considered  as  a  single  variable  type  in  this  research. 

Discrete  variables  are  common  in  engineering  design  problems.  For  example, 
in  the  design  of  aircraft,  the  number  of  engines  is  integer- valued  and  the  engine  type 
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Figure  1.2:  Graphical  Depiction  of  the  Objective  Space  (Non-Convex  Pareto  Front) 

(turboprop,  turbofan,  etc.)  and  engine  placement  (wing,  aft,  or  combination)  [126] 
are  categorical.  Other  discrete  variables  may  include  airfoil  type,  wing  conhguration 
(orientation  and  location  of  the  lifting  surfaces  of  the  aircraft  [8, 127]),  and  cabin 
layout  (number  of  aisles,  location  of  lavatories,  etc.  [127]).  The  class  of  optimization 
problems  that  may  contain  continuous,  discrete- numeric,  and  categorical  variables  is 
known  as  mixed  variable  programming  (MVP)  [9,138]. 

To  model  mixed  variables,  the  decision  space  is  partitioned  into  continuous  and 
discrete  domains,  and  respectively,  where  the  discrete  variables  may  include 
categorical  variables.  Without  loss  of  generality,  discrete  variables  can  be  mapped 
to  the  integers,  so  that  the  discrete  part  of  the  decision  space  can  be  represented  as  a 
subset  of  the  integers,  he.,  C  where  is  the  dimension  of  the  discrete  space. 
A  solution  X  E  Q  is  denoted  by  a;  =  (a;^,  a;'^),  where  x'^  G  K"’'"  and  x‘^  G  Z"'”^. 

1.1.4  Problem  Formulation.  The  inclusion  of  stochastic  and  multi-objective 
elements  to  the  classic  optimization  problem  formulation  yields  the  following  restated 
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problem: 


min  E[F{x)]  =  E[f{x)  +£a;(a^)]  (1.5a) 

subject  to 

gi{x)<0,  i  =  (1.5b) 

where  x  G  (M"''"  x  Z"‘*)  and  E  :  (R”'"  x  Z””^)  — >  R*^. 

1.2  Requirements  Analysis  and  Researeh  Outline 

Because  engineering  design  optimization  and  many  other  practical  optimization 
applications  are  generally  multi-objective  and  may  contain  both  stochastic  elements 
and  mixed  variables,  an  optimization  method  capable  of  handling  the  following  four 
characteristics  is  desired. 

1.  Stochastic.  Stochastic  optimization  algorithms  should  be  able  to  approxi¬ 
mate  the  Pareto  frontier  when  starting  from  an  arbitrary  point  for  problems  in 
which  function  evaluations  are  known  to  contain  measurement  error  or  must  be 
estimated,  e.g.,  via  simulation.  Further,  such  an  algorithm  should  be  provably 
convergent  to  Pareto  solutions,  thus  guaranteeing  that  a  representation  of  the 
frontier  in  a  region  of  interest  can  be  found.  However,  convergence  for  stochastic 
methods  is  usually  specihed  as  “almost  surely”  or  “with  probability  1”  [137]. 

2.  Multi-Objective.  Ideally,  algorithms  should  be  capable  of  hnding  a  rea¬ 
sonably  accurate  approximation  of  the  Pareto  frontier  in  cases  for  which  no 
preference  information  is  explicitly  known  or  even  exists.  However,  in  many 
engineering  design  problems,  some  information  about  desired  performance  goals 
(called  aspiration  levels),  as  well  as  minimum  acceptable  performance  levels 
(called  reservation  levels),  may  in  fact  exist,  and  can  be  used  to  determine  a 
region  of  interest  in  which  to  estimate  the  Pareto  front.  This  type  of  preference 
information  is  assumed  to  exist. 
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3.  General  Purpose.  An  algorithm  should  be  applicable  to  a  wide  range  of 
problems,  with  any  combination  of  variable  types.  An  algorithm  should  also  be 
indifferent  or  robust  to  the  source  of  the  function  evaluations;  he.,  the  algorithm 
is  able  to  treat  the  objective  and  constraint  functions  as  “black  boxes”. 

4.  Efficient.  To  be  practical  and  useful  for  real  design  problems,  an  algorithm 
should  perform  well  with  respect  to  the  number  of  function  evaluations  required. 
In  many  design  applications,  such  function  evaluations  are  obtained  via  costly 
simulation  runs  and  should  therefore  be  used  as  parsimoniously  as  possible. 

1.2.1  Problem  Statement.  There  exists  no  provably  convergent,  general- 
purpose  class  of  methods  for  solving  multi-objective,  stochastic  optimization  problems 
that  apply  to  the  mixed-variable  case  and  are  indifferent  to  the  source  of  function 
evaluations,  as  well  as  computationally  efficient  algorithmic  implementations  of  these 
methods. 

1.2.2  Research  Objectives.  The  purpose  of  this  research  is  to  develop  a 
provably  convergent,  general-purpose  class  of  methods  for  solving  multi-objective, 
stochastic  optimization  problems  that  apply  to  the  mixed-variable  case  and  are  in¬ 
different  to  the  source  of  function  evaluations.  Each  of  these  requirements  for  the 
solution  method  presents  specific  challenges  in  the  optimization  process.  In  response 
to  these  unique  challenges,  a  new  method  is  introduced,  which  extends  the  applica¬ 
bility  of  generalized  pattern  search  with  ranking  and  selection  for  linearly  constrained 
problems  [138]  and  mesh  adaptive  direct  search  (MADS)  for  nonlinearly  constrained 
problems  [6]  to  multi-objective  problems  through  the  use  of  interactive  specification 
of  aspiration/reservation  levels  [62],  scalarization  functions  [101],  and  multi-objective 
ranking  and  selection  methods  [85] .  Specific  goals  of  the  research  are: 

1.  Determine  an  appropriate  amalgamation  of  stochastic  and  multi-objective  meth¬ 
ods  to  be  extended  to  the  multi-objective,  stochastic  case.  Investigate  tradeoffs 
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among  the  desirable  characteristics  of  proven  convergence  properties,  simplicity, 
computational  efficiency,  and  processing  time. 

2.  Extend  previous  methods  to  address  the  mixed- variable,  multi-objective,  stochas¬ 
tic  class  of  problems.  Examine  and  prove  convergence  properties  of  the  new 
solution  methodology. 

3.  Develop  specihc  algorithmic  implementations  of  the  new  solution  methodology. 
Test  the  implementations  on  a  range  of  appropriate  test  problems  and  evaluate 
computational  efficiency.  The  following  sub-objectives  summarize  the  empirical 
results  pursued. 

(a)  Test  the  algorithm  on  a  set  of  test  problems  with  known  solutions.  In  order 
to  evaluate  the  algorithm’s  use  on  stochastic  problems,  the  test  problems 
are  modihed  to  introduce  random  noise  into  the  objective  function  evalua¬ 
tions.  Compare  the  solution  quality  and  the  computational  efficiency  to  the 
published  results  of  other  solution  methods.  Because  no  methods  currently 
exist  for  the  mixed-variable,  multi-objective,  stochastic  class  of  problems, 
compare  the  new  algorithm  to  those  used  for  deterministic,  multi-objective 
problems  and  infer  a  quality  evaluation  of  computational  speed  and  accu¬ 
racy. 

(b)  Test  the  algorithm  on  an  engineering  design  problem. 

4.  Study  appropriate  algorithm  termination  criteria.  Develop  heuristic  stopping 
criteria  based  on  differences  in  competing  responses  compared  to  variations  in 
the  responses  and  the  practically  required  tolerance  of  the  solution. 

1 . 3  Overview 

This  dissertation  is  organized  as  follows.  Chapter  II  reviews  several  existing 
solution  methods  for  both  stochastic  and  multi-objective  problems,  as  well  as  the 
few  that  exist  for  problems  with  both  characteristics.  Chapter  III  presents  specihc 
elements  of  the  proposed  methodology,  an  outline  of  the  proposed  methodology,  and 
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theoretical  convergence  resnlts.  Chapter  IV  provides  algorithmic  implementations 
of  the  solntion  methodology.  Chapter  V  presents  compntational  resnlts  of  these 
implementations  tested  against  mnlti-objective  test  problems  with  known  resnlts,  as 
well  as  an  engineering  design  optimization  problem.  Finally,  Chapter  VI  snmmarizes 
the  achievements  of  this  research  and  snggests  areas  for  farther  research. 
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II.  Literature  Review 


As  illustrated  in  Figure  2.1,  optimization  is  a  diverse  field  with  many  branches.  Prob¬ 
ably  the  most  well  known  and  longest  studied  problem  is  that  of  single  objective, 
deterministic  optimization.  Examples  of  traditional  optimization  techniques  include 
Dantzig’s  simplex  method  for  linear  programs  [162],  branch  and  bound  techniques 
and  decomposition  like  Sweeney-Murphy  decomposition  for  integer  programs  [142] 
and  Benders  decomposition  for  large-scale  problems  [21],  nonlinear  programming 
techniques  like  sequential  quadratic  programming  [73],  and  heuristic  methods  like 
evolutionary  algorithms.  However,  though  applicable  to  many  types  of  problems, 
other  methods  are  needed  for  problems  with  more  than  one  objective  or  for  systems 
having  stochastic  responses.  Thus,  the  disciplines  of  stochastic  optimization  and 
multi-objective  optimization  have  been  developed,  but,  for  the  most  part,  separately. 
The  following  sections  provide  an  overview  of  stochastic  and  multi-objective  optimiza¬ 
tion  methods,  as  well  as  methods  that  have  been  used  for  stochastic  multi-objective 
problems  with  continuous  variables. 

2.1  Methods  for  Stochastic  Optimization 

Stochastic  optimization  problems  are  those  in  which  the  objective  or  constraint 
functions  contain  a  random  element,  either  due  to  random  noise  in  the  function  eval¬ 
uations  or  random  algorithm  choices  (such  as  the  search  direction)  [134].  Stochastic 
optimization  problems  can  be  further  delineated  by  the  nature  of  solution  required. 
Parametric  optimization  is  used  to  find  values  for  a  set  of  parameters  to  optimize 
a  performance  measure,  e.g.^  cost,  reliability,  or  weight.  In  contrast,  control  opti¬ 
mization  seeks  to  find  a  set  of  actions  to  be  taken  in  different  states  a  system  can 
visit  in  order  to  optimize  a  performance  measure  [61].  A  classic  example  of  this 
type  of  problem  is  the  Markov  decision  problem,  which  is  typically  solved  using  dy¬ 
namic  programming.  Various  solution  techniques  for  both  types  of  problems  are 
shown  in  Figure  2.1.  (Note  that  when  multiple  objectives  can  be  aggregated  in  some 
way — see  Section  2.2.1 — single  objective  control  optimization  techniques  can  be  used 
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for  multi-objective  problems.)  Though  this  research  addresses  the  held  of  paramet¬ 
ric  optimizatiou,  a  brief  review  of  coutrol  optimizatiou  is  preseuted  for  the  sake  of 
completeuess  before  preseutiug  methods  of  parametric  optimizatiou. 

2.1.1  Solution  Methods  for  Control  Optimization. 

2. 1.1.1  Dynamic  Programming.  Developed  iu  the  late  1950’s  by  Bell- 
mau  aud  Howard  [20, 69] ,  dyuamic  programmiug  has  a  cousiderably  lower  computa- 
tioual  burdeu  thau  exhaustive  euumeratiou.  It  has  coutiuued  to  evolve  over  the  last 
half  ceutury  aud  is  uow  cousidered  to  be  the  maiu  pillar  of  stochastic  coutrol  opti¬ 
mizatiou.  Dyuamic  programmiug  assigus  a  so-called  value  function  to  the  process 
iu  each  state.  Theu  the  Bellman  equations  for  value  iteration — recursiou  equatious 
for  expected  reward  based  ou  the  probability  of  a  state  aud  the  value  associated  with 
that  state  [155] —  are  used  to  hud  the  optimal  set  of  actious  to  move  the  process  aloug 
the  best  path.  Simulatiou  cau  be  used  to  estimate  the  value  fuuctiou  at  each  state. 
If  the  coutrol  optimizatiou  problem  cousists  of  a  Markov  or  semi-Markov  [20, 69]  de- 
cisiou  process,  theu  dyuamic  programmiug  is  guarauteed  to  couverge  to  au  optimal 
solutiou. 


2. 1.1. 2  Reinforcement  Learning.  A  severe  drawback  to  dyuamic  pro¬ 
grammiug  is  that  the  trausitiou  probabilities  betweeu  states  must  be  kuowu  iu  order  to 
use  the  Bellmau  equatious,  which  may  be  difficult  to  obtaiu  for  real-world  problems. 
Overcomiug  this  problem  may  require  either  simulatiou  to  approximate  trausitiou 
probabilities  or  simply  usiug  a  method  that  does  uot  require  them.  Siuce  simula¬ 
tiou  is  somewhat  iuefficieut,  the  latter  choice  is  geuerally  preferred.  Reiuforcemeut 
learuiug  is  oue  such  method  [61]. 

Reiuforcemeut  learuiug  was  developed  iu  the  artihcial  iutelligeuce  commuuity 
aud  is  kuowu  as  a  “machiue  learuiug”  techuique.  Though  ouly  guarauteed  to  hud 
uear-optimal  solutious,  it  is  still  a  powerful  solutiou  method  because  it  is  based  ou 
the  Markov  decisiou  model,  lu  this  approach,  the  decisiou  maker  is  viewed  as  a 
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learning  agent,  who  receives  feedback  from  the  environment  based  on  decisions  the 
agent  makes.  Good  decisions  prodnce  reinforcement  by  increasing  the  probability 
the  agent  will  take  that  action,  bnt  bad  ones  resnlt  in  pnnishment  by  decreasing 
this  probability.  In  this  way,  the  agent  learns  by  trial  and  error  in  a  simnlated 
environment  [61].  Specihc  reinforcement  learning  methods  inclnde  adaptive  critic 
methods  [15],  Q-learning  [61,151],  and  temporal  difference  learning  [134]. 

2.1.2  Methods  for  Stochastic  Parametric  Optimization.  Consider  the  prob¬ 
lem  [61], 

/OO 

[x  —  5]^^'^p{x)dx, 

-OO 

where  p  is  a  probability  density  fnnction  (pdf)  of  the  random  variable  x,  and  the 
objective  fnnction  /  is  the  expected  valne  of  x.  Since  the  pdf  is  known  and  / 
is  known  and  has  closed  form,  this  problem  can  be  solved  with  standard  nonlinear 
programming  techniqnes  snch  as  steepest  descent,  gradient  projection  methods  for 
constrained  problems,  and  linearization  methods  like  Frank- Wolfe  [52,102,107].  If 
the  closed  form  of  the  objective  fnnction  is  not  known  or  has  no  closed  form,  or  if 
it  is  difficnlt  to  obtain  valnes  for  the  pdf,  other  methods  mnst  be  nsed.  Types  of 
solntion  methods  that  estimate  the  valne  of  the  objective  fnnction  nsing  simnlation  are 
aptly  called  simulation-based  optimization.  Simnlation,  thongh  not  an  optimization 
techniqne  per  se,  can  be  nsed  in  conjnnction  with  nnmerical  optimization  methods 
in  order  to  solve  difficnlt  parametric  optimization  problems  [61].  Sriver  provides  an 
excellent  overview  of  simnlation-based  optimization  techniqnes  [137].  These  (and 
additional)  methods  are  discnssed  in  Sections  2. 1.2. 1-2. 1.2. 3. 

2. 1.2.1  Response  Surface  Methodology.  Response  snrface  methodology 
(RSM),  developed  by  Box  et  al.  [26,27,104,105],  is  one  of  the  oldest  simnlation-based 
methods  of  parametric  optimization.  It  is  based  on  the  idea  that  an  approximate 
form  of  the  objective  fnnction  can  be  obtained  by  simnlating  the  system  at  a  finite 
nnmber  of  carefnlly  chosen  points  in  the  decision  space.  Traditionally,  RSM  has 
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used  polynomial  regression  over  the  sampled  points  to  find  an  approximate  form  of 
the  objective  function;  however,  alternate  forms  have  been  proposed  such  as  kriging, 
neuro-response  surfaces,  multivariate  adaptive  regression  splines  (MARS),  and  radial 
basis  functions  [61, 137, 138].  A  survey  of  such  methods  can  be  found  in  [75].  It  is 
decidedly  robust  and  often  works  well  where  other  methods  fail  [61];  however,  it  only 
converges  to  the  optimal  value  of  the  model,  not  that  of  the  original  problem. 

2. 1.2. 2  Continuous  Optimization.  Continuous  parametric  optimiza¬ 
tion  refers  simply  to  optimization  problems  with  continuous  variables.  This  class  of 
problems  has  two  general  types  of  solution  methods:  those  that  use  the  gradient  of 
the  objective  function  and  those  that  do  not  [61, 162]. 

Gradient-Based  Methods.  Stochastic  gradient-based  methods 
include  gradient  descent  and  simultaneous  perturbation.  First,  consider  gradient 
descent  (also  called  the  method  of  steepest  descent  or  finite  difference  stochastic  ap¬ 
proximation).  A  stochastic  version  of  one-dimensional  gradient  descent  for  finding 
the  root  of  a  one-dimensional  unconstrained  noisy  function  was  first  introduced  by 
Robbins  and  Monro  in  the  1950s  [123].  It  was  extended  shortly  thereafter  by  Keifer 
and  Wolfowitz  [76]  to  that  of  using  central  finite  differences  to  estimate  the  derivative, 
and  later  by  Blum  to  include  the  multivariate  case  [23] .  It  has  been  the  most  popular 
and  widely  used  method  for  simulation-based  optimization  [143]. 

At  each  iteration,  the  partial  derivatives  are  estimated,  the  new  search  direction 
is  chosen  as  the  negative  of  the  approximated  gradient,  and  the  algorithm  moves  a 
small  step  in  this  direction.  This  repeats  until  the  estimate  of  the  gradient  is  close 
to  zero  (within  a  tolerance). 

The  gradient  descent  method,  which  is  defined  by  the  iteration, 

=a;^-^V/(a;^),  =  1,2,...,  (2.1) 
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is  guaranteed  to  converge  to  a  stationary  point  if  the  function  /  is  convex,  continu¬ 
ously  differentiable,  and  the  step  size  ^  is  sufficiently  small  [61].  When  the  closed 
form  of  the  objective  function  is  not  known,  the  “sufficiently  small”  step  size  must 
be  found  by  experimentation.  Additionally,  2n  function  evaluations  are  required  per 
iteration  to  estimate  the  gradient  in  Equation  (2.1),  and  since  the  function  evalua¬ 
tions  are  generated  via  simulation,  multiple  replications  may  be  required  per  function 
evaluation.  Thus,  as  the  number  of  decision  variables  grows,  the  number  of  (possibly 
time-consuming)  simulation  runs  grows  dramatically  [61,134]. 

Spall  [134]  developed  simultaneous  perturbation,  a  method  in  which  the  gradient 
can  be  estimated  with  only  two  function  evaluations,  regardless  of  the  number  of 
variables.  The  objective  function  is  evaluated  at  the  current  iterate  and  at  another 
perturbed  point  a  small  distance  from  the  current  iterate  in  order  to  approximate 
partial  derivatives.  The  key  is  that  the  perturbed  point  is  shifted  (randomly)  in  all 
dimensions  at  once.  This  method  is  not  only  efficient,  but  also  has  been  shown  to 
converge  with  probability  1  to  a  stationary  point  of  the  objective  function  [61]. 

Gradient-Free  Methods.  The  second  type  of  continuous  opti¬ 
mization  solution  methods  do  not  require  gradients,  e.g.,  direct  search  methods  [137]. 
These  methods  evolved  from  early  efforts  to  optimize  stochastic  systems.  In  the 
1950s,  Box  proposed  a  method  for  improving  industrial  efficiency  known  as  evolu¬ 
tionary  operation  (EVOP)  [24]  in  which  estimation  of  regression  models  was  replaced 
by  a  relatively  simpler  method  of  inspecting  experimental  data  according  to  patterns 
in  the  experimental  design.  Later,  an  automated  version  of  EVOP  with  updated 
experimental  designs  (simplex  vice  factorial)  was  proposed  by  Spendley,  Hext,  and 
Himsworth  [135]. 

During  the  1960s,  several  more  direct  search  methods  were  developed  includ¬ 
ing  the  search  method  of  Hooke  and  Jeeves  [67]  and  the  Nelder-Mead  method  [109]. 
The  Nelder-Mead  method,  also  known  as  the  downhill  simplex  method  or  flexible 
polygon  search,  is  an  extension  of  the  method  developed  by  Spendley,  Hext,  and 
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Himsworth  [135].  These  two  methods  are  perhaps  the  most  frequently  used  direct 
search  methods  and  have  been  widely  used  on  real-world  problems  with  much  success. 
At  the  time,  they  were  generally  regarded  as  heuristics  because  they  were  thought  to 
not  possess  satisfactory  convergence  properties  in  high  dimensions  [61]  or  a  general 
convergence  theory  in  the  stochastic  setting  [137].  More  recently,  Torczon  proved 
convergence  for  a  more  general  class  of  pattern  search  methods,  which  includes  the 
direct  search  method  of  Hooke  and  Jeeves  but  not  the  Nelder-Mead  method  [146]. 
(Torczon  also  states  that  the  simplex  method  of  Spendley,  Hext,  and  Himsworth  can 
be  shown  to  converge  with  some  modihcation  to  the  original  algorithm,  but  only  under 
the  additional  assumption  that  /  is  convex  [146].)  Later,  McKinnon  [100]  showed  ex¬ 
amples  in  which  Nelder-Mead  converges  to  nonstationary  points.  These  examples  are 
for  well-behaved  functions  and  illustrate  the  limitations  on  the  convergence  properties 
of  Nelder-Mead  [100].  A  survey  of  pattern  search  methods  can  be  found  in  [80]. 

The  Nelder-Mead  method  starts  with  a  feasible  set  of  solutions  (dimension  nj-l) 
whose  vertices  dehne  a  “simplex”  or  “polygon” .  At  each  iteration,  the  worst  solution 
is  dropped  from  the  simplex  in  favor  of  a  good  one.  During  the  search,  the  geometry 
of  the  simplex  is  changed  by  expansion,  reflection,  contraction,  or  shrinking  operations 
based  on  the  relative  rank  of  the  point  added  at  each  iteration  [109].  Improvements 
and  extension  to  stochastic  problems  of  this  method  include  work  by  Tomick,  Arnold, 
and  Barton  to  determine  rules  for  estimating  the  number  of  replications  required  at 
each  simplex  design  point  [145];  a  variant  of  the  method  proposed  by  Barton  and  Ivy 
that  reduces  the  risk  of  false  convergence  [16];  and  work  by  Humphrey  and  Wilson 
that  addresses  sensitivity  to  starting  values,  premature  termination,  robustness,  and 
computational  efficiency  [71,72]. 

The  Hooke- Jeeves  pattern  search  method  involves  two  types  of  movement.  First, 
a  search  area  is  determined  via  exploratory  moves  from  the  starting  point.  A  pattern 
of  points  is  then  dehned  in  this  area  and  is  explored  by  changing  each  parameter  in 
turn.  The  pattern  is  then  shifted  to  a  new  region  via  the  exploratory  moves  [140]. 
This  method  has  been  used  in  conjunction  with  two-stage  ranking  and  selection  pro- 
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cedures  [111]  as  well  as  inside  the  optimization  modnles  of  the  GASP  and  SLAM 
simnlation  langnages  [113,114].  Little  research  was  done  on  direct  search  methods 
from  the  1960s  nntil  the  introdnction  and  convergence  analysis  of  generalized  pattern 
search  (GPS)  for  nnconstrained  optimization  problems  by  Torczon  in  the  1990s  [146]. 
GPS  is  discnssed  in  more  detail  in  Ghapter  111. 

2. 1.2. 3  Discrete  Optimization.  Discrete  optimization  problems  have 
gaps  in  the  domain  of  the  objective  fnnction,  so  derivatives  may  be  of  little  nse.  If 
the  solntion  space  is  snfhciently  small,  then  the  problem  can  be  solved  by  exhanstive 
search  of  the  solntion  space.  However,  as  problem  dimension  increases,  exhanstive 
search  qnickly  becomes  impractical.  Additionally,  even  when  the  solntion  space  is 
relatively  small,  a  more  efficient  method  for  comparing  and  ranking  the  best  solntions 
is  preferred.  [61] 


Ranking  and  Selection.  The  ranking  and  selection  techniqnes 
of  Goldsman  and  Nelson  [60]  or  Hochberg  and  Tamhane  [66]  are  examples  more 
efficient  methods.  Ranking  and  selection  methods  are  statistical  methods  for  selecting 
the  best  solntion  from  among  a  set  of  competing  candidates.  These  methods  are 
theoretically  sonnd,  well  tested,  and  are  typically  nsed  on  systems  with  np  to  twenty 
candidate  solntions,  althongh  recent  work  shows  that  they  can  be  nsed  on  mnch 
larger  problems  [61].  Other  examples  of  ranking  and  selection  algorithms  are  the 
Rinott  method  and  the  Kim-Nelson  method,  a  more  recent,  and  possibly  more  efficient 
algorithm  [61]. 


Meta-heuristics.  When  several  hnndred  or  thonsand  candidate 
solntions  exist,  ranking  and  selection  methods  cannot  be  nsed  directly.  They  can, 
however,  be  nsed  inside  of  a  meta-henristic  to  increase  the  power  of  the  ranking  and 
selection  method  (by  increasing  the  nnmber  of  alternatives  that  can  be  examined) 
and  to  make  the  resnlts  of  the  meta-henristic  more  reliable.  Meta-henristics  focns 
on  hnding  “good”  solntions  to  large-scale  (possibly  otherwise  intractable)  discrete 
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optimization  problems.  They  do  not  require  a  closed  form  of  the  objective  function, 
relying  only  on  value  of  the  function  at  various  points  which  can  be  estimated  via 
simulation.  While  they  are  not  guaranteed  to  produce  optimal  solutions,  they  often 
produce  good  solutions  in  a  reasonable  amount  of  time  and  are  included  in  some  com¬ 
mercial  simulation  software  packages.  Examples  of  meta-heuristics  include  simulated 
annealing,  genetic  algorithms,  tabu  search,  scatter  search,  and  the  learning  automata 
search  technique  [61, 134].  A  brief  description  of  each  method  follows. 

Simulated  annealing  is  an  adaptive  search  method  that  “exploits  the  analogy 
between  the  way  in  which  metal  cools  and  then  freezes  into  a  minimum  energy  crys¬ 
talline  structure  and  the  search  for  a  minimum  in  a  more  general  system”  [147].  In 
simulated  annealing,  search  directions  (or  points,  in  the  discrete  case)  are  randomly 
generated.  If  a  point  with  a  lower  function  value  is  found,  then  it  becomes  the  new 
iterate.  Otherwise,  it  may  still  randomly  accept  the  point  with  nonzero  probability, 
based  on  the  function  values  and  a  control  parameter.  While  it  may  seem  odd  to 
accept  a  worse  point,  doing  so  can  help  avoid  local  (non-global)  minima  [128,134,147]. 

Genetic  algorithms,  one  class  of  evolutionary  algorithms,  mirror  the  evolutionary 
processes  found  in  nature.  In  nature,  each  species  searches  for  benehcial  adaptation 
in  a  changing  environment.  As  species  evolve,  new  attributes  are  encoded  in  the 
chromosomes  of  its  members.  Though  this  information  can  change  by  mutation, 
most  change  occurs  due  to  the  combination  and  exchange  of  genetic  attributes  during 
breeding.  Analogously,  genetic  algorithms  code  complicated  non-biological  structures 
as  simple  bit  strings  and  improve  the  structures  by  simple  transformations  of  the 
strings.  They  are  unique  in  that  they  work  on  encoded  decision  variables  rather 
than  the  variables  themselves,  and  search  from  one  population  to  another  rather  than 
from  one  individual  to  another  [147].  Thus,  each  new  population  has  a  good  chance 
of  inheriting  the  best  characteristics  of  the  previous  one  [134].  Genetic  algorithms 
are  generally  effective,  robust,  computationally  simple  [128],  and  easy  to  parallelize 
[134, 147].  However,  they  are  often  very  expensive,  in  that  they  generally  require 
many  function  evaluations. 
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Tabu  search,  proposed  by  Glover  in  1977,  is  described  as  “a  meta-heuristic 
superimposed  on  another  heuristic”  [57].  The  crux  of  tabu  search  is  to  avoid  getting 
caught  in  cycles  or  local  optima  by  forbidding  or  penalizing  moves  which  steer  the 
solution  towards  regions  of  solution  space  already  visited.  As  with  other  meta¬ 
heuristics,  it  is  partially  motivated  by  an  observation  of  nature,  human  nature  in 
this  case.  Human  behavior  appears  to  operate  with  a  random  element  that  leads 
to  inconsistent  behavior  in  similar  situations.  It  is  this  inconsistent  behavior  or 
deviation  from  a  charted  course  that  proves  to  be  a  source  of  gain.  Tabu  search 
works  in  this  way  except  that  it  does  not  randomly  choose  new  courses.  Instead  it 
only  accepts  a  “poor”  course  if  it  is  the  only  way  to  avoid  a  path  previously  studied. 
In  this  way  it  can  potentially  escape  the  trap  of  local  minima  [128, 131]. 

Scatter  search  is  a  population-based  method  originally  proposed  by  Glover  in 
the  1960s.  Unlike  other  evolutionary  algorithms,  it  does  not  search  in  a  completely 
random  fashion.  It  instead  uses  strategies  for  combining  effective  solution  vectors 
through  the  use  of  an  adaptive  memory  similar  to  that  used  by  tabu  search  [18, 19, 
58,103]. 

The  learning  automata  search  technique  (LAST)  is  based  on  game  theory  and 
artihcial  intelligence.  Unlike  many  other  meta-heuristics,  LAST  is  not  a  local  search 
technique.  Instead,  it  jumps  around  the  solution  space  while  keeping  the  best  point. 
Initially,  it  jumps  randomly,  but  as  the  algorithm  progresses,  it  begins  to  “zero  in”  on 
its  perceived  best  solution  by  updating  the  probabilities  associated  with  the  random 
jump.  Initially,  the  algorithm  has  an  equal  probability  of  choosing  various  points 
to  evaluate.  At  each  iteration,  a  point  is  selected  based  on  probabilities  that  are 
stored  with  each  previously  evaluated  point,  and  evaluated.  If  the  selected  point  has 
a  comparatively  low  objective  function  value,  then  the  corresponding  probability  is 
rewarded  with  an  increase.  Otherwise,  it  is  punished  with  a  decrease.  In  this  way 
the  algorithm  “gets  smarter”  [61,164].  LAST  has  also  been  combined  with  a  genetic 
algorithm  to  speed  solution  convergence  properties  and  increase  its  chance  of  escaping 
local  optima  [70]. 
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2. 1.2. 4  Mixed  Variable  Solution  Spaces.  Often  in  real-world  problems 
the  solntion  space  is  not  neatly  defined  as  either  discrete  or  continnons,  bnt  rather  a 
combination  of  both  types,  i.e.,  it  is  a  problem  in  which  some  decision  variables  are 
discrete  and  some  are  continnons.  Historically,  methods  to  tackle  this  problem  have 
inclnded  1)  discretize  the  solntion  space  of  the  continnons  fnnctions  nsing  a  reasonably 
hne  grid,  2)  treat  the  discrete  variables  as  continnons  and  nse  continnons  optimization 
algorithms  with  ronnding,  or  3)  nse  branch-and-bonnd  techniqnes  branching  only  on 
the  discrete  variables.  Natnrally,  these  methods  have  drawbacks.  With  the  hrst 
option,  the  design  space  can  become  qnite  large  depending  on  the  size  of  the  grid, 
while  with  the  second,  treating  discrete  variables  as  continnons  is  known  to  have 
problems  in  practice,  often  leading  to  snb-optimal  or  infeasible  solntions.  Finally, 
branch-and-bonnd  techniqnes  reqnire  relaxation  of  integrality  constraints  [162]  and 
so  are  not  applicable  to  problems  involving  categorical  variables.  Solntion  methods 
for  this  type  of  problem  inclnde  pattern  search  methods  for  solving  deterministic 
mixed  variable  problems  introdnced  by  Andet  and  Dennis  [9].  Recent  work  by  Sriver 
et  al.  [138]  extended  the  work  of  Andet  and  Dennis  to  the  stochastic  case  nsing  a 
combination  of  pattern  search  and  ranking  and  selection.  This  method  (MGPS-RS) 
is  discnssed  fnrther  in  Chapter  111.  A  more  general  framework  for  treating  mixed 
variable  problems  (withont  the  specifying  the  approach  for  handling  the  continnons 
variables)  was  introdnced  by  Lncidi  et  al.  for  both  the  derivative-free  case  [95]  and 
the  case  in  which  derivatives  are  employed  [94]. 

2.2  Methods  for  Multi-objective  Optimization 

There  are  many  methods  available  to  solve  mnlti-objective  optimization  prob¬ 
lems.  These  methods  can  be  sorted  into  three  families:  a  priori  methods,  progressive 
methods,  and  a  posteriori  methods.  With  a  priori  methods,  the  decision  maker  de- 
hnes  the  trade-off  to  be  applied  (preferences)  before  rnnning  the  optimization  method. 
With  progressive  methods,  the  decision  maker  improves  the  trade-off  as  the  method 
progresses.  And  finally,  in  a  posteriori  methods,  the  decision  maker  chooses  a  solntion 
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after  examining  all  those  obtained  by  the  optimization  method.  Not  all  methods  fall 
neatly  into  one  family.  For  example,  a  priori  preferences  can  be  computed  at  random 
and  then  a  large  number  of  solutions  can  be  presented  to  the  decision  maker  to  choose 
the  hnal  solution. 

Most  multi-objective  optimization  methods  can  also  be  divided  into  hve  sets: 
scalar  methods,  interactive  methods,  fuzzy  methods,  meta-heuristic  methods,  and 
decision  aid  methods,  as  shown  in  Figure  2.1(a).  Again,  not  all  methods  can  be 
pigeon-holed  into  a  single  set.  For  example,  the  method  of  interactive  specihcation  of 
aspiration/reservation  levels,  discussed  extensively  in  Chapter  111,  is  a  scalar  method 
that  is  also  interactive.  Though  not  absolute,  the  sets  do  provide  a  useful  way  of 
categorizing  and  describing  most  methods.  A  description  and  examples  of  each  set 
of  methods  follows.  Additional  information  on  multi-objective  solution  methods  can 
be  found  in  [33]. 

2.2.1  Scalar  Methods.  Scalar  methods  attempt  to  transform  the  multi¬ 
objective  problem  into  one  with  a  single  objective,  so  that  standard  optimization 
techniques  can  be  applied.  Several  methods  are  presented. 

2. 2. 1.1  Weighted  Sum  of  the  Objective  Functions.  This  approach  is  the 
most  obvious  of  the  solution  methods  and  is  sometimes  called  the  “naive  approach” . 
It  simply  constructs  the  single  objective  as  a  weighted  sum  of  the  objective  functions. 
The  problem  can  be  solved  repeatedly  using  different  different  weights  to  generate  an 
approximation  of  the  Pareto  front.  This  was  the  first  solution  method  used  on  multi¬ 
objective  problems;  though  it  is  quite  efficient,  it  may  perform  poorly  on  problems 
with  nonconvex  Pareto  fronts.  If  the  Pareto  front  is  nonconvex,  not  all  solutions  on 
the  tradeoff  surface  can  be  found  [33, 37, 128] .  For  example,  the  Pareto  front  depicted 
in  Figure  1.2  is  an  extreme  case  where  only  the  two  endpoints  of  the  front  can  be 
found  with  this  type  of  method,  regardless  of  the  weights  [32] . 
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2. 2. 1.2  Keeney-Raiffa  Method.  The  Keeney- Raiffa  method  forms  a 
single  objective  fnnction  as  the  prodnct  of  the  objective  fnnctions.  The  new  objective 
function  is  called  the  Keeney-Raiffa  ntility  fnnction  and  is  nsed  in  mnlti-attribnte 
ntility  theory,  which  comes  from  the  decision  sciences  and  deals  with  the  properties 
of  a  ntility  fnnction  and  ways  to  create  it  [33]. 

2. 2. 1.3  Distance-to-a-Reference  Point  Method.  This  method  also 
transforms  the  objective  fnnctions  into  one  by  calcnlating  the  distance  to  a  refer¬ 
ence  point,  snch  as  the  ntopia  point,  nsing  an  appropriate  distance  metric  {e.g..,  the 
Enclidean  norm).  A  Pareto  optimal  point  is  then  fonnd  by  minimizing  the  distance 
to  the  reference  point.  The  snccess  of  this  method  is  heavily  dependent  on  the  choice 
of  a  reference  point,  bnt  it  is  often  able  to  hnd  solntions  hidden  in  concavities  of  the 
Pareto  front  [33]. 

2. 2. 1.4  Compromise  Method.  While  the  previons  three  methods  merge 
several  objective  fnnctions  into  a  single  one,  the  compromise  method  does  so  with 
additional  constraints.  This  is  done  by  preserving  the  highest  priority  objective 
fnnction  (as  specihed  by  the  decision  maker)  and  transforming  all  others  into  con¬ 
straints.  Thongh  the  compromise  method  is  known  to  consnme  large  amonnts  of 
compnting  time  and  the  programming  of  the  algorithm  can  be  difhcnlt  with  many 
objective  fnnctions,  the  relative  simplicity  of  its  eqnations  has  made  it  popnlar  [33]. 

2. 2. 1.5  Normal  Boundary  Intersection  Method.  The  normal  bonnd- 
ary  intersction  (NBl)  approach  of  Das  and  Dennis  [38]  prodnces  an  approximation 
of  the  Pareto  front  by  solving  a  series  of  single-objective  optimization  problems,  in 
which  an  additional  linear  eqnality  constraint  based  on  previonsly  determined  objec¬ 
tive  fnnction  valnes  is  added  [12].  This  resnlts  in  a  Pareto  solntion  that  lies  on  the 
intersection  of  the  constraint  and  the  Pareto  front.  It  is  based  on  the  observation 
that  the  intersection  of  the  bonndary  of  the  set  of  feasible  solntions  (or  image  of  the 
feasible  region)  with  a  normal  pointing  towards  the  origin  and  emanating  from  any 
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point  in  the  convex  hnll  of  individnal  minima  (convex  combination  of  the  objectives) 
is  a  point  on  the  portion  of  the  bonndary  which  contains  efficient  points.  NBI  has 
been  shown  to  provide  an  evenly  distribnted  set  of  points  on  the  Pareto  front  [37] . 

2. 2. 1.6  Hybrid  Methods.  Hybrid  methods  combine  several  methods 
into  new  ones.  The  best  known  hybrid  method  is  the  Corley  Method,  which  con¬ 
verts  some  objectives  to  constraints  and  then  applies  a  weighted-snm-of-the-objective- 
fnnctions  to  the  remaining  objective  fnnctions.  This  method  is  efficient  in  solving  var- 
ions  kinds  of  optimization  problems,  with  both  convex  and  nonconvex  Pareto  fronts, 
bnt  has  twice  the  nnmber  of  parameters  to  “tnne”  (e.g'.,  weights,  which  objectives  to 
convert  to  constraints,  and  in  what  order)  [33]. 

2. 2. 1.7  Goal  Attainment  Method.  The  goal  attainment  method  is  rel¬ 
atively  simple  and  is  able  to  determine  points  on  nonconvex  regions  of  the  Pareto 
front  (see  Fignre  1.1).  The  methods  consists  of  choosing  an  initial  vector  of  decision- 
maker-specihed  ideal  fnnction  valnes,  choosing  a  search  direction  (essentially  weights 
or  relative  importance  of  the  objectives),  and  minimizing  a  scalar  coefficient  which 
represents  the  gap  between  the  two  vectors.  The  main  advantage  of  this  method 
is  its  efficiency  with  respect  to  nonconvex  Pareto  fronts;  however,  some  shapes  of 
Pareto  fronts  exist  for  which  the  snccess  of  the  method  depends  greatly  on  the  choice 
of  the  search  direction  [33] .  An  example  is  the  method  of  interactive  specihcation  of 
aspiration/reservation  levels,  mentioned  in  Section  2. 2. 2. 6,  which  is  an  interactive/- 
goal  attainment  method  that  nses  aspiration  levels  (decision-maker-specihed  ideal  or 
“hoped  for” )  as  the  reference  point.  The  search  direction  is  determined  by  the  ray  be¬ 
tween  the  aspiration  level  and  reservation  level  (decision-maker-specihed  “worst  case 
scenario”)  [96].  This  method  is  discnssed  fnrther  in  Section  3.1.4. 

2. 2. 1.8  Goal  Programming  Method.  This  method  is  similar  to  the 
goal  attainment  method,  except  that  the  transformed  problem  has  eqnality  instead 
of  ineqnality  constraints.  Becanse  it  is  so  similar,  this  method  has  the  same  relative 
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advantages  and  disadvantages  as  the  goal  attainment  method  [33]. 

2. 2. 1.9  Lexicographic  Method.  As  the  name  implies,  this  method  takes 
the  objectives  and  orders  them  by  preference.  The  mono-objective  problems  are 
then  solved  one  at  a  time  while  converting  already  solved  objectives  into  constraints. 
Thongh  intnitive  and  easily  solved,  the  main  drawback  of  this  method  is  that  an 
appropriate  seqnence  of  the  objective  fnnctions  mnst  be  determined  withont  inadver¬ 
tently  exclnding  promising  regions  of  the  space  prematnrely  [33, 128]. 

2.2.2  Interactive  Methods.  Interactive  methods  belong  to  the  progressive 
methods  family  and  thns  allow  the  decision  maker  to  tnne  preferences  with  regard  to 
tradeoffs  as  the  methods  progress.  These  methods  hnd  one  and  only  one  solntion. 
Several  interactive  methods  are  presented  below. 

2.2.2. 1  Surrogate- Worth  Tradeoff  Method.  This  method  adds  an  in¬ 
teractive  process  to  the  compromise  method  (see  Section  2.2. 1.4)  to  progress  towards 
a  solntion.  Since  the  decision  maker  has  inpnt  thronghont  the  process,  a  solntion 
fonnd  via  this  method  has  a  good  chance  of  being  satisfactory.  However,  the  choice  of 
opinion  valnes  can  greatly  affect  the  solntion  qnality,  and  becanse  it  relies  on  gradient 
information,  this  method  shonld  not  be  applied  to  problems  for  which  the  objective 
fnnctions  are  not  differentiable  [33,106]. 

2.2. 2. 2  Fandel  Method.  The  goal  of  the  Fandel  method  is  to  assist  the 
decision  maker  in  choosing  the  weights  nsed  in  the  weighted-snm-of-the-objective- 
fnnctions  method.  Each  objective  fnnction  is  hrst  optimized  separately  to  generate 
the  ntopia  point.  This  point  and  the  feasible  region  are  then  nsed  to  calcnlate  a 
domain  of  variation  (region  of  the  objective  space  bonnded  by  the  minimnm  and 
maximnm  valnes  of  each  objective).  An  attempt  is  then  made  to  hnd  the  hyperplane 
parallel  to  the  one  passing  throngh  the  extremities  of  the  domain  of  variation  and 
tangent  to  the  feasible  space.  The  tangent  point  provides  a  vector  that  corresponds 
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to  weights  that  are  used  to  determine  a  solution.  The  decision  maker  either  accepts 
this  solution  or  the  search  space  is  reduced.  The  main  drawback  of  the  Fandel  method 
is  that  parts  of  the  tradeoff  surface  cannot  be  approximated,  and  it  only  hnds  a  single 
solution.  Additionally,  it  assumes  that  the  solution  that  satishes  the  decision  maker 
is  the  one  closest  to  the  ideal  point,  which  may  not  be  the  case  [33]. 

2. 2. 2. 3  Step  Method.  This  method  is  similar  to  the  Fandel  method  in 
that  information  about  the  preferences  of  the  decision  maker  allows  the  search  space 
to  be  restricted  at  each  iteration.  The  decision  maker  indicates  an  upper  bound 
for  acceptable  values  of  the  objective  functions  and  then  weights  are  calculated  such 
that  objectives  with  a  larger  range  of  possible  values  are  given  larger  weights  in  the 
weighted-sum-of-the-objective-functions.  The  problem  is  optimized  and  the  new 
solution  is  evaluated  by  the  decision  maker.  If  it  is  not  acceptable,  the  upper  bound 
of  the  solution  space  can  be  adjusted  and  the  weights  recalculated.  Compared  with 
the  Fandel  method,  the  weights  have  less  importance  because  interaction  with  the 
decision  maker  is  done  via  selection  of  the  desired  upper  bound  and  not  the  weights 
directly  [33]. 

2. 2. 2. 4  Jahn  Method.  In  contrast  to  the  Fandel  and  STEP  methods, 
the  Jahn  method  depends  on  both  the  decision  maker  preferences  and  gradient  infor¬ 
mation.  It  begins  by  choosing  a  starting  point  and  then  stepping  through  the  solution 
space  using  a  gradient-based  method  and  decision  maker  preferences  toward  a  Pareto 
optimal  point.  Since  gradient  information  is  required,  the  objective  functions  must 
be  differentiable  [33]. 

2.2.2. 5  Geoffrion  Method.  This  method  is  similar  to  the  Jahn  method 
and  is  based  on  the  Frank- Wolfe  algorithm.  The  method  begins  with  a  weighted- 
sums-of-the-objective-function  step.  The  resulting  point  is  then  used  to  linearize  the 
objective  function  and  compute  a  direction  that  reduces  the  objective  while  main¬ 
taining  feasibility  [52].  In  the  next  step  instead  of  performing  a  line  search  to  find 
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the  optimal  step  along  that  direction  (as  in  the  Frank- Wolfe  method),  the  step  size 
is  determined  by  the  decision  maker.  The  new  solntion  is  either  accepted  by  the 
decision  maker  or  the  process  is  repeated  [33]. 

2. 2. 2. 6  Interactive  Specification  of  Aspiration/Reservation  Levels  and 
Achievement  Scalarization  Functions.  As  mentioned  in  Sections  2.2  and  2. 2. 1.7, 
this  method  is  based  on  idea  that  a  decision  maker  has  a  desired  level  for  each  ob¬ 
jective,  as  well  as  a  level  beyond  which  the  solntion  is  not  acceptable.  These  valnes, 
called  aspiration  levels  and  reservation  levels,  respectively,  are  nsed  in  a  method  sim¬ 
ilar  to  distance-to-a-point  or  compromise  methods  (where  the  aspiration  level  is  the 
reference  point  and  both  are  nsed  to  determine  the  search  direction)  in  order  to  hnd  a 
single  solntion.  The  decision  maker  can  then  either  accept  the  solntion  or  revise  the 
aspiration  and  reservation  levels  [96].  This  method  is  discnssed  fnrther  in  Chapter 
111. 


2.2.3  Fuzzy  Methods.  Thongh  real  world  phenomenon  are  rarely  binary  in 
natnre,  for  a  long  time  the  only  tool  for  describing  them  was  binary  logic,  written 
in  terms  of  true  or  false.  In  response,  L.A.  Zadeh  developed  a  new  logic  based  on 
fuzzy  sets,  which  not  only  offers  a  way  to  model  nncertainty  and  imprecision,  bnt  also 
acconnts  for  progressive  transition  between  states,  e.g.  shades  of  gray  between  black 
and  white.  This  is  accomplished  with  a  membership  function,  which  relates  each 
decision  variable  to  a  continnons  variable  varying  between  zero  and  one  and  indicates 
the  degree  to  which  the  variable  belongs  to  the  set  represented  by  the  original  logical 
decision  variable  and  allowing  a  progressive  transition  between  true  and  false.  Two 
fuzzy  logic  methods  are  briefly  presented  [33,128]. 

2.2.3. 1  Sakawa  Method.  This  method  nses  fnzzy  logic  at  all  levels: 
the  parameters  of  the  problem,  constraints,  and  solntion  set.  Typical  of  fnzzy  logic 
methods,  the  solntions  have  a  membership  fnnction  that  is  correlated  with  the  orig¬ 
inal  objective  fnnction  at  a  level  set  by  the  decision  maker.  This  method  allows 
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the  decision  maker  to  specify  the  number  of  solutions  desired,  which  is  useful  when 
the  decision  maker  is  more  demanding  [33].  Sakawa  also  combined  fuzzy  logic  with 
evolutionary  methods  in  order  to  handle  multi-objective  binary  programming  prob¬ 
lems  [129].  Despite  its  interesting  advantages,  the  Sakawa  method  is  cumbersome 
and  uses  fuzzy  sets  and  mathematical  rules  that  are  difficult  to  apply  [33]. 

2. 2. 3. 2  Reardon  Method.  In  contrast,  the  Reardon  method  is  a  rela¬ 
tively  simplihed  fuzzy  logic  method.  Its  membership  function  is  simply  shaped  and 
easily  obtained  compared  to  the  Sakawa  method.  It  has  been  shown  to  give  good 
results  on  two  test  problems,  Schaffer  F2  test  problem  and  the  simplified  Born-Mayer 
test  problem  [33]  and  that,  in  contrast  to  genetic  algorithms  like  the  niched  Pareto 
approach  (see  Section  2. 2. 4. 4),  the  efficiency  of  the  algorithm  is  independent  of  the 
number  of  objectives  [120, 121]. 

2.2.4  Multi- objective  Methods  Using  Meta- Heuristics.  Meta- heuristics  are 

general  optimization  methods  which  need  some  transformation  before  being  applied 
to  the  solution  of  a  particular  problem  and  are  often  analogous  to  real  life  concepts 
[128].  They  are  usually  applied  to  “hard”  optimization  problems.  As  discussed 
in  Section  2. 1.2. 3,  meta-heuristics  are  not  guaranteed  to  produce  optimal  solutions, 
but  they  often  produce  good  solutions  in  a  reasonable  amount  of  time  [33].  Four 
meta-heuristics  are  presented:  simulated  annealing,  tabu  search,  scatter  search  and 
genetic  (evolutionary)  algorithms. 

2.2.4. 1  Simulated  Annealing.  Simulated  annealing,  discussed  pre¬ 
viously  for  stochastic  optimization  in  Section  2. 1.2.3,  has  also  been  used  in  multi¬ 
objective  optimization  [141].  Examples  of  multi-objective  simulated  annealing  meth¬ 
ods  include  Multiple  Objective  Simulated  Annealing  (MOSA)  [33,47,93],  Fuzzy  Pareto 
Simulated  Annealing  (FPSA)  [133],  and  Pareto  Archived  Simulated  Annealing  (PASA) 
[33].  MOSA  uses  simulated  annealing  to  search  the  tradeoff  surface,  considering  ob¬ 
jectives  separately  to  determine  the  probability  of  accepting  a  bad  solution  but  using 
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a  weighted-sum-of-the-obj  active- functions  approach  to  calculate  the  actual  solution 
htness  (he.  the  value  of  the  weighted  objective  function).  The  method  behaves  well 
because  at  a  high  temperature,  the  simulated  annealing  spreads  the  solution  popu¬ 
lation  over  the  whole  tradeoff  surface  [33].  If  the  temperature  is  “sufficiently  high” 
at  the  end  of  the  simulation,  MOSA  can  find  solutions  unreachable  by  the  weighted- 
sum-of-the-obj active- functions  method  [147];  however,  what  constitutes  “sufficiently 
high”  is  completely  problem-dependent  is  not  generally  known  a  priori.  In  FPSA, 
simulated  annealing  is  combined  with  fuzzy  logic  (see  Section  2.2.3)  in  order  to  find  an 
approximation  of  the  Pareto  set  [65].  PASA  uses  an  aggregated  function  of  objective 
functions,  coupled  with  a  system  that  stores  nondominated  solutions  separately  from 
the  rest  of  the  iterates  and  uses  them  to  assess  solution  quality  [33]. 

2. 2-4 .2  Tabu  Search.  Tabu  search,  discussed  previously  for  stochastic 
optimization  (Section  2. 1.2.3),  has  been  successfully  used  in  multi-objective  searches 
for  continuous  problems  [33,55].  Implementations  include  M-Objective  Tabu  Search 
[28],  MOTS  [47,64],  and  New  Multi-objective  Tabu  Search  (NMTS)  [118].  Tabu 
search  has  also  been  used  in  conjunction  with  scatter  search  algorithms,  e.g.,  Scatter 
Tabu  Search  Procedure  for  Multiobjective  Optimization  (SSPMO)  [103]. 

2. 2. 4. 3  Scatter  Search.  Scatter  search,  discussed  previously  for  stochas¬ 
tic  optimization  (Section  2. 1.2. 3),  has  been  successfully  used  in  multi-objective  searches 
[18,58,97,98].  Implementations  include  the  Scatter  Tabu  Search  Procedure  for  Mul¬ 
tiobjective  Optimization  (SSPMO)  [103]  and  Multi-Objective  Scatter  Search  (MOSS) 

|19]. 


2. 2. 4. 4  Genetic/Evolutionary  Algorithms.  Genetic  algorithms,  a  type 
of  evolutionary  algorithm,  were  also  discussed  as  a  stochastic  optimization  method  in 
Section  2. 1.2. 3.  They  are  useful  for  multi-objective  problems  because  they  deal  simul¬ 
taneously  with  a  set  of  possible  solutions  (the  population)  which  allows  the  algorithm 
to  find  several  elements  of  the  Pareto  optimal  set  in  a  single  run.  Unlike  many  tra- 
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ditional  math  programming  techniques  {e.g.,  aggregation  methods),  they  can  handle 
discontinuous  and  concave  Pareto  fronts.  A  multi-objective  genetic  algorithm  differs 
from  single-objective  ones  in  how  they  evaluate  solutions.  Unlike  single-objective  ge¬ 
netic  algorithms,  multi-objective  ones  cannot  directly  use  the  objective  function  value 
as  a  measure  of  solution  fitness,  so  some  combination  of  the  objectives  must  be  used. 
Many  methods  have  been  proposed,  and  multi-objective  evolutionary  algorithms  can 
be  classified  by  how  they  determine  if  a  solution  is  efficient:  aggregating  functions, 
population-based  approaches,  and  Pareto-based  approaches  [32,134]. 

Aggregating  Functions.  As  discussed  in  Section  2.2.1,  aggregating 
functions  are  probably  the  most  straightforward  of  approaches,  but  certain  problems 
exist.  Linear  aggregating  functions,  a  weighted  sum  for  example,  cannot  generate 
nonconvex  portions  of  the  Pareto  front  no  matter  what  weights  are  chosen.  How¬ 
ever,  nonlinear  functions — e.g.,  achievement  scalarization  functions  (see  Section  3.1.4 
and  [96]),  the  Keeney-Raiffa  method  (see  Section  2. 2. 1.2),  and  the  goal  attainment 
method  (see  Section  2. 2. 1.7) — do  not  suffer  from  this  limitation  and  have  been  used 
successfully  in  multi-objective  combinatorial  optimization  [32]. 

Population- Based  Approaches.  Population-based  approaches  use 
the  population  of  the  evolutionary  algorithm  to  diversify  the  search  and  potentially 
find  multiple  Pareto  optimal  points  at  each  iteration,  but  the  idea  of  dominance  is 
not  directly  incorporated  into  the  selection  process.  An  example  of  this  type  of 
approach  is  the  Vector  Evaluated  Genetic  Algorithm  (VEGA)  developed  by  Schaf¬ 
fer  [132].  VEGA  consists  of  a  simple  genetic  algorithm  with  a  modified  selection 
mechanism.  At  each  generation,  k  sub-populations  of  size  N/k,  where  N  is  the  total 
population  size,  are  generated  based  on  proportional  selection  according  to  each  objec¬ 
tive  function,  in  turn  [132].  These  sub-populations  are  then  shuffled  together  to  form 
one  population  to  which  crossover  and  mutation  operators  are  applied.  VEGA  has 
several  problems,  the  most  serious  of  which  is  that  the  selection  method  it  uses  does 
not  accurately  incorporate  Pareto  dominance.  If,  for  example,  a  good  compromise 
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(Pareto  optimal)  solution  is  generated,  but  it  is  not  the  best  solution  for  any  given 
objective  function,  it  will  be  discarded.  Some  methods  have  been  developed  to  help 
alleviate  this,  but  the  underlying  problem  remains  [32,47].  Other  population-based 
approaches  include  the  Random  Sampling  Evolutionary  Algorithm  (RAND)  [165]  and 
Hajela  and  Lin’s  genetic  algorithm  (HLGA)  [63,165]. 

Pareto-Based  Approaches.  In  contrast  to  population-based  ap¬ 
proaches,  this  type  of  approach  incorporates  the  concept  of  Pareto  optimality  di¬ 
rectly  into  the  selection  method  of  the  evolutionary  algorithm  by  using  dominance 
rather  than  an  aggregated  objective  function  to  measure  htness.  Thus,  the  Pareto- 
based  approaches  do  not  suffer  the  problems  associated  with  VEGA.  Many  Pareto- 
based  approaches  have  been  developed  [32],  including  Goldberg’s  Pareto  Ranking 
[59],  Multi- Objective  Genetic  Algorithm  (MOGA)  [47,51],  the  Nondominated  Sort¬ 
ing  Genetic  Algorithm  (NGSA  and  NGSA-11)  [42,136,165],  Niched  Pareto  Genetic 
Algorithm  (NPGA)  [68, 165],  Strength  Pareto  Evolutionary  Algorithm  (SPEA  and 
SPEA2)  [47,165],  Multi-Objective  Messy  Genetic  Algorithm  (MOMGA,  MOMGA- 
11,  and  MOMGA-lla)  [40],  Multi-Objective  Hierarchical  Bayesian  Optimization  Al¬ 
gorithm  (hBOA)  [32],  Pareto  Archived  Evolutionary  Strategy  (PAES)  [47,79,112], 
Micro-Genetic  Algorithm  for  Multi-Objective  Optimization  [31],  and  the  General 
Multi- Objective  Program  (GENMOP)  [78,125]. 

2.2.5  Decision  Aid  Methods.  Decision  aid  methods  differ  from  the  others 
in  several  ways.  Instead  of  filtering  the  elements  of  the  solution  set  and  keeping 
those  elements  that  could  be  compared  to  themselves  {e.g.,  lexicographic  optimality), 
an  order  relation  is  set  up  among  elements  of  the  solution  set  and  thus  a  set  of 
solutions  (with  a  partial  order  relation,  i.e.,  a  relation  that  is  reflective,  asymmetric, 
and  transitive  [153]),  or  a  single  solution  (with  a  total  order  relation,  i.e.,  a  relation 
that  is  reflective,  asymmetric,  transitive,  and  totaT  [154,160]),  is  obtained.  Also, 

binary  relation  R  over  a  set  X  is  total  if  it  holds  for  all  a  and  6  in  X  that  a  is  related  to  b  or 
b  is  related  to  a  (or  both),  i.e.,  Va,  b  G  X,  aRb  V  bRa  [161]. 
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these  methods  only  work  on  discrete  sets  of  points.  The  decision  aid  methods  choose 
or  sort  actions  with  respect  to  a  set  of  criteria.  Each  method  prodnces  its  own 
definition  of  the  valne  of  a  criterion  [33].  Decision  aid  methods  inclnde  ELECTRE 
methods  (I,  IS,  II,  III,  IV,  and  TRI)  [43,56,130]  and  PROMETHEE  (I  and  II)  [130]. 

2.3  Methods  for  Stochastic  Multi- objective  Optimization 

Generally,  there  are  two  types  of  stochastic  optimization  problems:  those  with 
random  noise  in  the  fnnction  evalnations,  and  those  with  random  algorithm  choices 
(snch  as  the  search  direction)  [134].  Most  stochastic  mnlti-objective  problems  solved 
to  date  are  of  the  second  type.  Examples  inclnde  methods  that  nse  evolntionary 
algorithms  (see  Section  2. 2. 4. 4).  Relatively  few  mnlti-objective  solntion  methods  of 
the  hrst  type  have  been  developed.  Three  methods  proposed  by  Baba  and  Morimoto 
[13]  are  discnssed  in  Sections  2.3.1,  2.3.2,  2.3.3.  A  more  recent  method  developed  by 
Andet  et  al.  [12]  is  discnssed  in  Section  2.3.4.  This  method  has  not  yet  been  applied 
to  the  stochastic  case. 

2.3.1  Learning  Automata.  Learning  antomata  is  a  reinforcement  (or  feed¬ 
back)  learning  scheme  where  actions  by  the  antomaton  prodnce  resnlts  for  which 
either  reward  or  pnnishment  ensne.  The  feedback  then  changes  the  probability  of 
choosing  that  action  snch  that  rewards  increase  the  probability  of  selecting  that  ac¬ 
tion  again,  and  pnnishment  decreases  that  probability.  This  process  is  repeated  nntil 
the  antomaton  “learns”  what  actions  are  optimal  [61, 134].  In  the  mnlti-objective 
case.  Baba  and  Morimoto  show  that  an  appropriately  chosen  learning  scheme  ensnres 
convergence  to  a  “reasonable  solntion”  for  a  hnite  nnmber  of  candidate  solntions  [13] . 

2.3.2  Random  Optimization  Method.  Random  search  methods — typically 
nsed  for  single  objective  problems — are  based  on  searching  the  problem  domain  in  a 
random  manner  in  order  to  optimize  an  objective.  These  methods  are  the  simplest 
of  stochastic  optimization,  bnt  can  be  effective  in  certain  cases  [134].  Baba  and 
Morimoto  showed  that  their  random  optimization  algorithm  is  gnaranteed  to  converge 
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almost  surely  to  a  Pareto-optimal  solution  but  only  under  strict  assumptions  on  the 
decision  space,  solution  space,  and  error.  They  suggest  further  study  to  hud  a  less 
restrictive  approach  [13]. 

2.3.3  Stochastic  Approximation  Method.  Baba  and  Morimoto  [13,14]  pro¬ 

pose  a  stochastic  quasigradient  method  to  solve  the  stochastic  multi-objective  opti¬ 
mization  problem.  Subgradients  are  used  for  convex  functions  that  are  not  necessarily 
differentiable  at  all  points.  For  the  convex  function  /  :  ^  R,  the  set  of  subgradients 

at  Xn  in  the  open  interval  is  Gn  =  {c  G  R  :  f{x)  —  f{xn)  >  c(x  —  Xn),  Va;  G  11}  [159]. 
A  stochastic  quasigradient  is  a  stochastic  estimate  of  the  subgradient  of  the  function 
at  the  point  Xn,  i.e.,  given  a  set  of  points  {xq  . . .  a;„},  the  stochastic  subgradient  is 
a  random  vector  such  that  £'(^„|a;o  . . .  Xn)  =  OnF^ixu)  +  &n  where  the  vector  Fx(xn) 
is  the  subgradient  of  the  function  F(x)  at  the  point  x  =  x^,  bn  is  a  7i- valued  ran¬ 
dom  variable,  Ti  is  a  separable  Hilbert  space,  and  a„  is  a  positive  real-valued  random 
variable.  [48,49,116].  Under  assumptions  of  continuity,  compactness,  and  bounded 
error.  Baba  and  Morimoto  show  that  the  algorithm  extended  to  the  multi-objective 
case  converges  with  probability  one  to  the  global  solution  [14] . 

2.3.4  Bi-Objective  Mesh  Adaptive  Direct  Search  (BiMADS).  Audet  et  al. 
have  extended  the  Mesh  Adaptive  Direct  Search  (MADS)  (discussed  in  Section  3.1.3) 
framework  to  apply  to  the  bi-objective  case.  Their  algorithm,  BiMADS,  is  similar  to 
reference  point  methods  (see  Section  2.2. 1.3),  but  also  takes  advantage  of  the  ordering 
property  of  the  Pareto  front  inherent  only  in  bi-objective  problems  [12,22].  BiMADS 
is  shown  to  produce  solutions  satisfying  some  first  order  necessary  conditions  for 
optimality  based  on  the  Clarke  calculus  [30],  and  like  the  NBI  method  (see  Section 
2. 2. 1.5),  attempts  to  generate  evenly  distributed  points  on  the  Pareto  front  [12].  Thus 
far  BiMADS  has  been  developed  for  deterministic  problems;  however,  the  extensions 
of  MADS  to  the  mixed  variable  and  stochastic  case  discussed  in  Chapter  III  would 
also  apply  to  BiMADS. 
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2-4  More  Generalized  Stochastic  Multi- objective  Optimization  Method 

As  seen  throughout  this  chapter,  many  solution  methods  exist  for  both  stochas¬ 
tic  optimization  and  multi-objective  optimization  when  only  one  of  the  two  charac¬ 
teristics  exists  in  the  problem.  As  discussed  in  Section  2.3,  four  methods  exist  for 
problems  having  both  characteristics  and  because  of  necessary  assumptions  on  the 
decision  and  objective  spaces,  such  methods  are  applicable  only  to  a  small  number  of 
problems.  None  of  these  methods  is  applicable  to  stochastic,  multi-objective,  mixed- 
variable  optimization  problems  as  considered  in  this  research.  Thus,  a  more  generally 
applicable  method  is  required  and  is  presented  in  Chapter  111. 
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III.  Research  Methodology 


3.1  Components  for  a  New  Algorithm 

In  this  chapter,  an  algorithm  is  presented  for  solving  stochastic  mixed  variable 
optimization  problems  with  mnltiples  objectives.  This  algorithm  relies  on  a  nnm- 
ber  of  important  featnres  that  are  described  here  in  greater  detail.  Particnlarly, 
consider  the  following  observations.  Generalized  pattern  search  with  ranking  and 
selection  (MGPS-RS)  has  been  snccessfnlly  developed  for  stochastic,  linearly  con¬ 
strained  mixed-variable  problems  [137, 138]  and  has  been  applied  to  a  mnlti-echelon 
repair  system  [144].  Mesh  adaptive  direct  search  (MADS)  has  been  developed  for  de¬ 
terministic,  continnons,  nonlinearly  constrained  problems  [11].  However,  MGPS-RS 
and  MADS  in  their  cnrrent  forms  apply  only  to  single-objective  problems.  Alter¬ 
natively,  interactive  techniqnes  nsing  aspiration/reservation  levels  and  scalarization 
functions  have  been  nsed  snccessfnlly  to  hnd  Pareto  optimal  solntions  to  determin¬ 
istic  mnlti-objective  problems  [62].  Finally,  a  mnlti-objective  ranking  and  selection 
techniqne  by  Lee  et  al.  [85],  called  multi-objective  optimal  computing  budget  allocation 
(MOGBA),  has  been  applied  to  selecting  the  non-dominated  set  of  inventory  policies 
for  aircraft  maintenance,  a  discrete  variable  problem  [86]. 

The  rest  of  this  section  gives  a  more  detailed  description  of  each  of  the  methods 
mentioned  thns  far.  A  new  algorithm,  called  Stochastic  Multi- Objective  Mesh  Adaptive 
Direct  Search  (SMOMADS),  is  then  proposed  in  Section  3.2,  which  combines  elements 
of  these  methods  to  handle  the  class  of  mixed  variable,  stochastic,  mnlti-objective 
optimization  problems.  A  convergence  analysis  of  the  new  algorithm  then  follows  in 
Section  3.3. 

3.1.1  Ranking  and  Selection.  Problems  having  stochastic  responses  can  pose 
a  particnlar  challenge  for  optimization  algorithms  that  rely  solely  on  comparisons  of 
trial  points.  Noise  in  the  response  can  lead  to  errors  in  iterate  selection  if  the  observed 
responses  are  comparatively  different  than  the  trne  fnnction  valnes.  Methods  for 
selecting  a  “best”  among  several  trial  points  mnst  acconnt  for  variation  to  provide 
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some  statistical  assurance  of  correct  selection.  Achieving  high  conhdence  in  the 
correctness  of  selection  requires  additional  replications  (he.,  function  evaluations) 
at  each  trial  point.  Ranking  and  selection  (R&S)  considers  multiple  candidates 
simultaneously  at  a  reasonable  cost.  This  is  possible  because  R&S  detects  a  relative 
order  of  the  candidates  instead  of  generating  precise  estimates  [139] .  The  discussion 
of  R&S  that  follows,  including  the  notation,  mirrors  that  of  [139]. 

Let  Xk  denote  the  k-th  element  of  a  sequence  of  random  vectors  and  Xk  denote 
a  realization  of  Xk-  For  a  finite  set  of  candidate  points  C  =  {Ti,y2;  •  •  •  with 

nc  >  2,  let  fq  =  f  (Yq)  =  E  [F  (Yq,  ■)]  denote  the  true  mean  of  the  response  function  F 
at  Yq  for  each  q  =  1,2, ...  ,nc-  These  response  means  can  be  ordered  (from  minimum 
to  maximum)  as 

/[l],/[2],  •  •  •  ,/[nc]-  (3-1) 

Denote  by  G  C  the  candidate  from  C  with  the  g-th  lowest  true  objective  function 
value. 

Given  some  <5  >  0,  called  the  indifference  zone  parameter,  no  distinction  is  made 
between  two  candidate  points  whose  true  means  satisfy  /p]  —  /[i]  <  <5.  In  such  cases, 
the  method  is  said  to  be  indifferent  in  choosing  either  candidate  as  the  best.  The 
probability  of  correct  selection  (CS)  is  dehned  as 

P{CS}  =  P  {select  Y[i]|/[g]  -  /[i]  >  5]  q  =  1,2,..  .,nc}  >  1  -  a,  (3.2) 

where  a  G  (0, 1)  is  the  signihcance  level.  Because  random  sampling  guarantees  that 

P  (C'S'}  =  — ,  the  signihcance  level  must  satisfy  0  <  a  <  1 - . 

nc  nc 

Because  the  true  objective  function  values  are  not  available,  it  is  necessary 
to  work  with  the  sample  means  of  F.  For  each  q  =  1,2, ...  ,nc,  let  Sq  be  the 
total  number  of  replications  at  Yq  and  let  =  {F  {Yqs,  kFgs)}^li  be  the  set  of 

simulated  responses,  where  are  the  replications  at  candidate  point  Y^,  and 

Wqs  are  realizations  of  the  random  noise.  For  each  q  =  1,2, ...  ,nci  the  sample  mean 
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Procedure  RS(C,  a,  <5) 

Inputs:  C  =  {Fi,  F2,  •  •  • ,  ^  (0, 1) ,  <5  >  0. 

Step  1  For  each  Yq  E  C,  use  an  appropriate  statistical  technique  to  de¬ 
termine  the  number  of  samples  Sq  required  to  meet  the  probability  of 
correct  selection  guarantee  in  Equation  (3.2),  as  a  function  of  a,  S  and 
respoinse  variation  of  Yq. 

Step  2  For  each  q  =  l,2,...,nc',  obtain  replicated  responses  Fqs,  s  = 

1,  2, . . . ,  Sq,  and  compute  the  sample  mean  Fq,  according  to  Equation 

(3.3). 

Return:  yp]  =  ai'g  (F|i|). 

Figure  3.1:  A  generic  R&S  Procedure  (see  Figure  1  in  [139]) 

Fq  is  given  by 

_  2  . 

Fq  =  ~'^^Fqs.  (3.3) 

*9  S=1 

The  sample  means  can  be  ordered  and  indexed  in  the  same  manner  as  in  Equation 
(3.1),  and  let  l^qj  G  C  denote  the  candidate  with  the  g-th  lowest  estimated  objective 
function  value  as  determined  by  the  R&S  procedure.  The  candidate  corresponding 
to  the  minimum  mean  response,  =  arg  (F[i]),  is  chosen  as  the  best  point.  A 
generic  R&S  procedure  is  shown  in  Figure  3.1.  The  algorithmic  choices  made  in  Step 
1  dehne  different  R&S  methods.  See  [110,115,122],  for  examples. 

3.1.2  MGPS-RS.  Pattern  search  algorithms  are  defined  through  a  finite  set 
of  directions  used  at  each  iteration.  The  direction  set  and  a  step  length  parameter 
are  used  to  construct  a  discrete  set  of  points,  or  mesh,  around  the  current  iterate. 
The  mesh  at  iteration  k  is  defined  to  be 

Mfc  =  [j  {x  +  A^Dz  :  ^  e  ,  (3.4) 

where  Sk  is  the  set  of  points  where  the  objective  function  /  has  been  evaluated  by 
the  start  of  iteration  k  (Sp  is  the  set  of  initial  points).  A™  is  called  the  mesh  size 
parameter,  and  D  is  a  positive  spanning  set  (see  [39]);  i.e.,  a  set  of  directions  that 
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positively  spans  Additional  restrictions  on  D  are  that  each  direction  d  ^  D, 

j  =  1,2, .. .  ,nD,  mnst  be  the  prodnct  of  some  hxed  nonsingnlar  generating  matrix 
G  G  and  an  integer  vector  zj  G  Z”"  [146].  For  bonnd  and  linearly  constrained 

problems,  the  directions  in  mnst  D  be  snlRciently  rich  to  ensnre  that  polling  directions 
can  be  chosen  that  conform  to  the  geometry  of  the  constraint  bonndaries,  and  that 
these  directions  be  nsed  inhnitely  many  times  [2],  A  hnite  set  of  trial  points  called 
the  poll  set  is  then  chosen  from  the  mesh,  evalnated,  and  compared  to  the  incnmbent 
solntion.  If  improvement  is  fonnd,  the  incnmbent  is  replaced  and  the  mesh  is  retained 
or  coarsened  via  the  mesh  size  parameter  A™.  If  not,  the  mesh  is  rehned  and  a  new 
set  of  trial  points  is  selected. 

Initially  developed  by  Torczon  [146] ,  GPS  was  extended  by  Lewis  and  Torczon  to 
inclnde  bonnd  [87]  and  linear  [88]  constraints.  GPS  was  farther  extended  by  Lewis 
and  Torczon  [89]  to  inclnde  nonlinear  constraints  nsing  an  angmented  Lagrangian 
approach  by  Gonn,  Gonld,  and  Toint  [34],  and  by  Andet  and  Dennis  [10],  via  a 
hlter  [50] .  Dissatisfaction  with  the  convergence  theory  of  hlter-GPS  led  to  the  creation 
of  MADS  [11]  (see  section  3.1.3).  GPS  for  bonnd  constrained  mixed  variable  problems 
(or  mixed  variable  pattern  search  -  MVPS)  was  introdnced  separately  by  Andet  and 
Dennis  [9] ,  and  extended  to  linearly  constrained  by  problems  by  Abramson  [2] .  It  was 
then  extended  fnrther  to  general  nonlinear ly  constrained  MVP  problems,  also  via  a 
hlter,  by  Abramson,  Andet  and  Dennis  [2,6]. 

The  GPS  framework,  in  conjnnction  with  ranking  and  selection,  was  nsed  by 
Sriver  [137, 139]  to  handle  stochastic  MVP  problems.  In  this  case,  the  poll  set  at 
each  iteration  is  given  by  Pk{xk)  IJ  A/"(a;fc),  where  J\f(xk)  is  a  nser-dehned  set  of  discrete 
neighbors  aronnd  Xk  and 


Pk{x)  =  {x  +  A^{d,0):deDl},  (3.5) 

where  {d,  0)  denotes  that  the  variables  have  been  partitioned  into  continnons  and  dis¬ 
crete  variables,  with  the  discrete  variables  remaining  nnchanged.  The  set  of  discrete 
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neighbors  is  defined  by  a  set-valned  fnnction  Af  :  ft  ^  2^,  where  2^  denotes  the  power 
set  of  fl.  The  notation  y  G  Af(x)  means  that  the  point  y  is  a  discrete  neighbor  of  x. 
By  convention,  x  G  Af(x)  for  each  a;  G  12,  and  it  is  assnmed  that  Af(x)  is  finite  [137]. 
A  generic  indifference-zone  ranking  and  selection  procednre  RS(Pk(xk)  IJ  A/"(a;fc),  a,  <5), 
with  indifference-zone  parameter  <5  and  significance  level  a,  is  nsed  to  select  among 
points  in  the  poll  set  for  improved  solntions,  i.e.,  (5-near-best  mean.  The  rnles  for 
npdating  the  mesh  size  are  as  follows  [11].  Given  a  fixed  rational  nnmber  r  >  1 
and  two  integers  m”  <  —  1  and  >  0,  the  mesh  size  parameter  A™  is  npdated 
according  to  the  rnle, 

A- 1  =  r-'=A-,  (3.6) 

where 

{{0, 1, ... ,  if  an  improved  mesh  point  is  fonnd 

(3.7) 

{m”,  -k  1, . . . ,  — 1},  otherwise. 

If  no  improvement  can  be  fonnd,  an  extended  poll  step  is  condncted  to  search  amongst 
any  discrete  neighbor  y  G  M{xk)  that  satisfies  f{y)  <  f(xk)  +  Ck,  where  is  called 
the  extended  polled  trigger.  Each  neighbor  meeting  this  criteria,  in  tnrn,  becomes 
the  poll  center,  and  the  extended  poll  continnes  nntil  either  a  better  point  than  the 
best  iterate  is  fonnd,  or  else  they  are  all  worse  than  the  extended  poll  center  [9, 137] . 
Sriver  showed  that  this  algorithm  generates  an  iteration  snbseqnence  with  almost 
snre  convergence  to  a  stationary  point  “appropriately  defined”  in  the  mixed-variable 
domain  [138].  The  MGPS-RS  Algorithm  is  shown  in  Fignre  3.2. 

3.1.3  Mesh  Adaptive  Direct  Search.  Mesh  Adaptive  Direct  Search  (MADS) 
is  a  class  of  algorithms  developed  by  Andet  and  Dennis  for  minimization  of  nonsmooth 
fnnctions  of  the  type  /  :  — >  IRlJ{-|-oo}  nnder  general  constraints  a;  G  D  C  R"-  and 
D  7^  0.  The  feasible  region  Q  may  be  defined  by  blackbox  constraints,  e.g.,  compnter 
code  that  retnrns  a  yes/no  answer  to  whether  or  not  a  trial  point  is  feasible  [11]. 
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A  General  MGPS-RS  Algorithm 

•  INITIALIZATION:  Let  Xq  G  G,  AJf  >  0,  ^  >  0,  ao  G  (0, 1)  and  (5o  >  0. 

Set  the  iteration  and  R&S  connters  /c  =  0  and  r  =  0  respectively. 

•  POLL  STEP:  Set  extended  poll  trigger  Use  R&S  procednre 

RS(Pk(Xk){jAf(Xk),ar,Sr)  to  retnrn  the  estimated  best  solntion  Y. 
Update  ar+i  <  ccr,  (5r+i  <  <5^,  and  r  =  r  +  l.  \i  Y  ^  X^,  the  step 
is  successful,  npdate  X^+i  =  Y,  >  A™  according  to  Eqnations 

(3.6)  and  (3.7),  and  k  =  k  +  1  and  retnrn  to  POLL  STEP.  Otherwise, 
proceed  to  EXTENDED  POLL  STEP. 

•  EXTENDED  POLL  STEP:  For  each  discrete  neighbor  Y  G  J\f{Xk) 
that  satishes  the  extended  poll  trigger  condition  F{Y)  <  F{Xk)  + 
set  j  =  1  and  Y^  =  Y  and  do  the  following. 

—  Use  R&S  procednre  RS(Pk(Y/!),ar,Sr)  to  retnrn  the  estimated 
best  solntion  Y.  Update  ar+i  <  ar,  (5r+i  <  (5r,  and  r  =  r  +  1. 

If  T  7^  yI,  set  yI'^^  =  Y  and  j  =  j  +  1  and  repeat  this  step. 
Otherwise,  set  Zk  =  Yf.  and  go  to  the  next  step. 

—  Use  R&S  procednre  RS{Xk[_}  Zk),oir,5r)  to  retnrn  the  estimated 
best  solntion  Y .  Update  <5^+1  <  <5,,,  and  r  =  r  +  1. 

If  T  =  Zk,  the  step  is  snccessfnl,  npdate  Xk+i  =  Y,  A^^  >  Xff 
according  to  Eqnations  (3.6)  and  (3.7),  and  k  =  k  +  1  and  retnrn 
to  the  POLL  STEP.  Otherwise,  repeat  the  EXTENDED  POLL 
STEP  for  another  discrete  neighbor  that  satishes  the  extended 
poll  trigger  condition.  If  no  snch  discrete  neighbors  remain  in 
N{Xk),  set  Xk+i  =  Xk,  <  A™  according  to  Eqnations  (3.6) 
and  (3.7),  and  k  =  k  +  1  and  retnrn  to  the  POLL  STEP. 

Fignre  3.2:  The  Mixed-variable  GPS  Ranking  and  Selection  (MGPS-RS)  Algorithm 
[138] 
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MADS  is  similar  to  MGPS-RS  in  the  generation  of  the  mesh  and  poll  sets  (see 
Eqnations  (3.4)  and  (3.5)  in  Section  3.1.2)  as  well  as  in  the  rnles  for  npdating  the  mesh 
(see  Eqnations  (3.6)  and  (3.7)).  However,  the  key  difference  is  that  in  MADS  [11]  a 
poll  size  parameter  is  introdnced,  which  satisfies  A™  <  A^  for  all  k  snch  that 

limA^  =  0  hmA|  =  0  for  any  infinite  snbset  of  indices  in  K.  (3.8) 

fceic  fcgic 

The  poll  size  parameter  controls  the  magnitnde  of  the  distance  between  the  incnmbent 
solntion,  and  the  set  of  directions  nsed  to  define  the  poll  set  are  constrncted  in  a 
different  manner.  In  GPS,  only  one  valne  A^  =  A^  =  A^  is  nsed,  and  a  set  of 
positive  spanning  directions  C  D  is  chosen  at  each  iteration.  In  the  poll  step  of 
MADS,  both  of  these  restrictions  are  relaxed,  in  that  neither  A^  =  A^  nor  C  D 
generally  hold.  The  poll  step  evalnates  points  lying  in  a  frame  (analogons  to  the 
poll  set  in  GPS),  constrncted  differently  than  GPS.  The  MADS  frame  is  defined  as 
follows  (see  [11]). 

Definition  3.1.1.  At  each  iteration  k,  the  MADS  frame  is  defined  to  he  the  set 

Pk  =  {xu  +  A^d  :  d  e  Dk]  C  Mk  (3.9) 

where  is  a  positive  spanning  set  such  that  0  ^  and  for  each  d  E  the  following 
conditions  must  be  met  [11]: 

1.  d  can  he  written  as  a  nonnegative  integer  combination  of  the  directions  in  D  : 
d  =  Du  for  some  vector  u  G  that  may  depend  on  the  iteration  number  k, 

2.  the  distance  from  the  frame  center  Xk  to  a  frame  point  Xk  +  A'^d  G  Pk  is  bounded 
above  by  a  constant  times  the  poll  size  parameter: 

Aff  II  d  ||<  A|max{|[  d'  |[:  d'  G  D} , 

3.  limits  of  the  normalized  sets  Dk  =  :  d  E  Dk\  are  positive  spanning  sets. 
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A  General  MADS  Algorithm 

•  INITIALIZATION:  Let  Xq  G  D,  <  Aq,  D,  G,  r,  w~ ,  and  satisfy 
the  reqnirements  of  a  MADS  frame  set  given  in  Dehnition  3.1.1. 

Set  the  iteration  connter  /c  ^  0. 

•  SEARCH  AND  POLL  STEP:  Perform  the  SEARCH  and  possibly  the 
POLL  steps  (or  part  of  them)  nntil  an  improved  mesh  point  Xk+i  is 
fonnd  on  the  mesh  (where  is  dehned  as  for  GPS  in  Eqnation 
(3.4)  in  Section  3.1.2). 

—  OPTIONAL  SEARCH:  Evalnate  /q  on  a  hnite  snbset  of  trial 
points  on  the  mesh  M^. 

—  LOCAL  POLL:  Evalnate  /q  on  the  frame  Pk  (where  Pk  is  as  given 
in  Eqnation  (3.9)  in  Section  3.1.2). 

•  PARAMETER  UPDATE:  Update  according  to  Eqnations  (3.6) 

and  (3.7)  and  ^k+i  Eqnation  (3.8)  is  satished. 

Set  /c  /c  +  1  and  go  back  to  the  SEARCH  AND  POLL  step. 

Fignre  3.3:  A  General  MADS  Algorithm  [11] 

The  general  MADS  algorithm,  as  developed  by  Andet  and  Dennis,  is  shown  in  Fignre 
3.3.  This  algorithm  is  extended  to  the  mixed  variable  case  in  [5]  and  in  Section  3.3.1. 
Extension  to  the  stochastic  case  is  accomplished  via  ranking  and  selection  in  a  manner 
similar  to  MGPS-RS.  The  extended  algorithm,  mixed  variable  mesh  adaptive  direct 
search  with  ranking  and  selection  (MVMADS-RS),  is  shown  in  Fignre  3.4. 

3.1.4  Interactive  Specification  of  Aspiration/Reservation  Levels  and  Scalariza- 
tion  Functions.  As  shown  in  Fignre  3.5(a),  points  on  the  Pareto  front  can  be 
fonnd  by  varying  the  relative  importance,  i.e.,  trade-off  coefficients  or  weights,  of  the 
distance  to  a  given  point.  Using  the  ntopia  point  U,  any  point  between  points  D 
and  E  can  be  fonnd.  (Points  D  and  E  correspond  to  Pareto  optimal  solntions  fonnd 
by  nsing  the  ntopia  point  in  one  dimension  and  the  nadir  point  in  the  other,  i.e.,  the 
best  possible  solntion  for  one  objective,  bnt  not  for  the  other.)  By  nsing  aspiration 
point  A  and  varying  the  weights  or  slope  of  the  ray  emanating  from  it,  points  be¬ 
tween  B  and  C  can  be  fonnd.  There  are  many  methods  for  determining  which  ray 
to  nse  [92] .  This  particnlar  method  nses  the  reservation  point  R  as  the  second  point 
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A  General  MVMADS-RS  Algorithm 

•  INITIALIZATION:  Let  Xq  G  O,  Af  >  A™  >  0,  ^  >  0,  ao  e  (0, 1) 
and  (5o  >  0.  Set  the  iteration  and  R&S  connters  /c  =  0  and  r  =  0 
respectively. 

•  POLL  STEP:  Set  extended  poll  trigger  Use  R&S  procednre 

RS(Pk(Xk)[jM(Xk),ar,Sr)  to  retnrn  the  estimated  best  solntion  Y. 
Update  ttr+i  <  ttr,  (5r+i  <  (5r,  and  r  =  r  +  l.  li  Y  ^  X^,  the  step  is 
successful,  npdate  X^+i  =  Y,  >  A|  according  to  Eqnations  (3.6) 
and  (3.7),  A^j^  >  A™  so  that  Eqnation  (3.8)  is  satished,  and  k  =  k  +  1 
and  retnrn  to  POLL  STEP.  Otherwise,  proceed  to  EXTENDED  POLL 
STEP. 

•  EXTENDED  POLL  STEP:  For  each  discrete  neighbor  Y  G  J\f(Xk) 
that  satishes  the  extended  poll  trigger  condition  F(Y)  <  F[Xjf)  + 
set  j  =  1  and  Yj!  =  Y  and  do  the  following. 

—  Use  R&S  procednre  RS(Pk(Yf!),ar,Sr)  to  retnrn  the  estimated 
best  solntion  Y.  Update  <5^+1  <  <5^,  and  r  =  r  +  1. 

If  T  7^  yI,  set  yI^^  =  Y  and  j  =  j  +  1  and  repeat  this  step. 
Otherwise,  set  Zk  =  Y/!  and  go  to  the  next  step. 

—  Use  R&S  procednre  RS{Xk[J  Zk),ar,Sr)  to  retnrn  the  estimated 
best  solntion  Y.  Update  a^+i  <  Q^r,  (5r+i  <  (5r,  and  r  =  r  +  1.  If 
Y  =  Zk,  the  step  is  snccessfnl,  npdate  X^+i  =  Y,  A^^^  >  A^  ac¬ 
cording  to  Eqnations  (3.6)  and  (3.7),  A^;^  >  A^  so  that  Eqnation 
(3.8)  is  satished,  and  k  =  k  +  1  and  retnrn  to  the  POLL  STEP. 
Otherwise,  repeat  the  EXTENDED  POLL  STEP  for  another  dis¬ 
crete  neighbor  that  satishes  the  extended  poll  trigger  condition. 
If  no  snch  discrete  neighbors  remain  in  jY{Xk),  set  X^+i  =  A^, 
<  Al  according  to  Eqnations  (3.6)  and  (3.7),  A^j^  <  A™  so 
that  Eqnation  (3.8)  is  satished,  and  k  =  k  +  1  and  retnrn  to  the 
POLL  STEP. 


Fignre  3.4:  The  Mixed-variable  MADS  with  Ranking  and  Selection  (MVMADS-RS) 
Algorithm 
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(a)  Pareto  solutions  corresponding  to  different 
component  achievement  functions  (Figure  3  in 
[96]) 


Ui 


(b)  Component  Achievement  Functions  for  Min¬ 
imized  Criteria  (Figure  4  in  [96]) 


Figure  3.5:  Graphical  illustrations  of  functions  used  for  analysis  of  Pareto  optimal 
solutions 


in  determining  the  direction  of  the  ray  [96]. 

This  technique  is  based  on  the  assumption  that  the  decision  maker  has  an  idea 
of  what  is  desired  for  each  objective,  as  well  as  what  minimum,  or  maximum,  values 
are  acceptable.  These  values  are  referred  to  as  the  aspiration  and  reservation  values, 
respectively;  he.,  points  A  and  R  discussed  previously  and  shown  in  Figure  3.5(a). 
These  values  are  then  used  inside  of  an  achievement  scalarization  function  of  the  form 
shown  graphically  in  Figure  3.5(b)  and  given  by 
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The  function  Ui  is  called  a  component  achievement  function,  he.,  strictly 
function  of  the  objective  function  vector  components  qi  =  fi{x),  i  =  1,2, . . 
by 
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where  g*  and  are  the  aspiration  and  reservation  levels,  respectively,  for  objective 
Wi  =  1/  (^q.  —  q^ ,  and  ttj,  Pi,  i  =  1,2, ...  ,J,  are  given  parameters  that  are  set  in 
snch  a  way  that  Ui  takes  the  valnes  at  the  ntopia,  aspiration,  reservation,  and  nadir 
points,  given  by 

respectively,  and  P  and  fj  are  given  positive  constants,  typically  set  to  0.1  and  10,  re¬ 
spectively  [96].  The  maximization  of  Eqnation  (3.10)  provides  proper  Pareto  optimal 
solntions  nearest  the  aspiration  level  (see  point  K  in  Fignre  3.5(a)). 

3.1.5  Multi- Objective  Ranking  and  Selection.  Lee  et  al.  [83]  propose  a  perfor¬ 
mance  index  to  measnre  the  degree  that  a  point  is  dominated  in  the  Pareto  sense  when 
the  objective  fnnction  evalnations  are  snbject  to  noise.  They  then  nse  this  measnre 
in  the  Mnlti-objective  Optimal  Compnting  Bndget  Allocation  (MOCBA)  algorithm 
to  determine  a  Pareto  optimal  set  for  mnlti-objective  simnlation-based  optimization 
problems  [86].  Given  a  set  of  designs  i  =  1,2, ...  ,n,  evalnated  by  J  performance 
measnres  fJ,ij,  j  =  1,2, . . . ,  J,  throngh  some  simnlation  (or  other  estimation  method), 
fj,ij  is  a  random  variable  whose  Bayesian  posterior  distribntion  can  be  derived  based 
on  its  prior  distribntion  and  the  simnlation  ontpnt  [83] .  Dehne  a  performance  index 

n 

'^i=  n  p  ii^k  <  tip] ,  (3.12) 

k=l,k^i 

where  P  {fXk  -<  t^p  represents  the  probability  that  design  k  dominates  design  i.  Under 
the  assnmption  that  the  performance  measnres  are  independent  and  follow  continnons 
distribntions  snch  that 
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ijji  measures  the  probability  that  design  i  is  not  dominated  by  any  other  designs 
and  can  be  used  inside  of  a  ranking  and  selection  framework  to  hnd  the  set  of  non- 
dominated  points  rather  than  a  single  best  point.  When  approximating  a  Pareto 
set,  two  types  of  error  exist:  Type  1  (ei),  the  probability  that  at  least  one  design  left 
out  of  the  efficient  set  is  non-dominated,  and  Type  11  (62),  the  probability  that  at 
least  one  design  in  the  efficient  set  is  dominated.  As  these  errors  both  approach  zero, 
the  experimental  Pareto  set  approaches  the  true^  Pareto  set.  It  is  shown  that  these 
errors  are  bounded  by  the  approximated  errors  aci  given  by 

ei  <  aci  = 

i€Sp 

and  ae2  given  by 

62  <  062  =  ^  (1  -  '4’i)  , 

i€Sp 

respectively,  and  that  the  observed  Pareto  set  determined  by  MOCBA  approaches 
the  true  Pareto  set  asymptotically  with  probability  1  [84, 85] . 

In  each  iteration  of  MGPS-RS  or  MADS,  only  a  hnite  number  of  points  are 
evaluated  by  the  ranking  and  selection  procedure.  Thus,  this  type  of  ranking  and 
selection  method  could  be  substituted  for  the  single  objective  ranking  and  selection 
method  inside  the  MGPS-RS  or  MADS  algorithms  to  develop  multi-objective  versions. 
Alternatively,  the  performance  measure  could  be  substituted  for  the  objective  func¬ 
tion  value  inside  of  a  traditional  ranking  and  selection  procedure.  In  this  case,  the 
convergence  results  from  [84, 85]  would  have  to  be  reverihed. 

^In  practice,  the  complete  true  Pareto  set  is  not  found.  In  this  context,  true  Pareto  set  refers 
to  an  approximation  containing  only  efficient  points  in  which  no  efficient  points  were  considered  but 
then  not  selected  as  such. 
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3.2  Stochastic  Multi- Objective  Mesh  Adaptive  Direct  Search 
(SMOMADS) 

A  new  two-stage  algorithm  which  uses  the  methods  of  Section  3.1  to  numer¬ 
ically  solve  mixed  variable  stochastic  multi-objective  optimization  problems  is  now 
introduced.  In  the  hrst  stage,  a  convex  combination  of  objectives,  via  scalarization 
functions  and  aspiration/reservation  levels  of  the  decision  maker,  is  used  to  deter¬ 
mine  an  approximation  of  the  Pareto  front  in  a  region  of  interest.  For  each  pair  of 
aspiration  and  reservation  levels,  MGPS-RS  or  MVMADS-RS  can  be  used  to  gen¬ 
erate  a  single  Pareto  point.  (Extension  of  MADS  to  the  mixed- variable  case — MV- 
MADS — is  discussed  in  Sections  3.2. 1.4  and  3.3.1  and  extension  to  the  stochastic 
case — MVMADS-RS — is  discussed  in  Section  3.3.2.)  However,  since  convexity  of  the 
actual  Pareto  frontier  is  not  assumed,  some  regions  in  the  Pareto  frontier  may  not 
contain  any  of  the  generated  points.  In  the  second  stage,  the  single-objective  rank¬ 
ing  and  selection  routine  inside  of  MGPS-RS  is  replaced  with  MOGBA,  so  that  the 
discrete  points  in  the  mesh  can  be  evaluated  with  respect  to  multiple  objectives.^  A 
graphical  representation  is  shown  in  Figure  3.6  and  descriptions  of  each  step  follow. 

3.2.1  Stage  1.  Stage  1  of  the  method  consists  of  aspiration  and  reservation 
level  analysis  in  order  to  represent  the  original  problem  as  one  with  a  single  objective 
that  can  then  be  solved  using  MGPS-RS  or  MVMADS-RS. 

3. 2. 1.1  Aspiration  and  Reservation  Level  Analysis.  As  discussed  in 
Section  3.1.4,  the  multiple  objectives  are  combined  into  a  single  objective  problem  of 
the  form  shown  in  Equation  (3.10).  Each  subproblem,  or  choice  of  a  specihc  pairing  of 
aspiration  and  reservation  levels,  produces  an  approximate  Pareto  point.  There  are 
many  ways  to  produce  such  aspiration/reservation  level  pairings.  (In  this  algorithm, 
the  problem  of  determining  these  pairings  is  referred  to  as  the  master  problem.)  His- 

^The  theoretical  convergence  properties  developed  in  Section  3.3.3  do  not  depend  on  stage  2  of 
the  algorithm.  It  was  added  to  better  determine  non-convex  Pareto  frontiers.  As  the  algorithm  was 
implemented  and  tested,  the  good  performance  of  Stage  1  made  Stage  2  unnecessary  for  the  problems 
tested.  Since  it  may  still  be  helpful  for  many  problems,  Stage  2  remains  part  of  the  algorithm. 
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Figure  3.6:  Stochastic  Multi-Objective  Mesh  Adaptive  Direct  Search  (SMOMADS) 


torically,  in  interactive  specihcation  of  aspiration  and  reservation  levels,  a  decision 
maker  was  actively  involved  in  choosing  these  points  [62],  However,  if  this  inter¬ 
action  is  not  possible  or  if  the  decision  maker  has  only  specihed  a  range  of  values 
for  aspiration  and  reservation  levels,  some  other  method  must  be  used.  In  the  case 
where  a  range  of  values  has  been  specihed,  the  problem  is  that  of  determining  an 
approximation  to  the  Pareto  frontier  within  a  region  of  interest.  Such  a  problem  is 
similar  to  that  of  approximating  a  response  surface  with  aspiration  and  reservation 
levels  as  the  decision  variables.  Thus,  experimental  design  methods  apply. 


3.2. 1.2  Example.  Given  a  problem  of  the  form 


min(/i(a:),/2(x)) 

where  D  represents  some  feasible  region  and  that  the  (as  yet  unknown)  Pareto  front 
is  as  shown  by  the  curve  in  Figure  3.7.  The  utopia  and  nadir  points  are  determined 
to  be  (1,3)  and  (6,8),  and  depicted  by  points  U  and  R,  respectively,  in  the  hgure. 
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Figure  3.7:  Graphical  Example  of  Aspiration/Reservation  Level  Analysis 


After  discussion  with  the  decision  maker,  it  is  decided  that  the  algorithm  should 
examine  aspiration  levels  between  1  and  3  for  objective  1  and  between  3  and  6  for 
objective  2  and  examine  reservation  levels  between  4  and  7  for  objective  1  and  between 
7  and  10  for  objective  2.  (This  is  depicted  by  the  rectangles  in  Figure  3.7.)  A  central 
composite  design  with  4  factors  (aspiration  and  reservation  levels  for  2  objectives)  is 
constructed,  resulting  in  36  test  points,  be.,  36  pairs  of  aspiration  and  reservation 
levels.  One  of  these  pairings  is  shown  by  the  points  A  and  R  in  Figure  3.7.  This 
pair  is  used  to  determine  the  component  achievement  functions  using  Equation  (3.11), 
which  become 


ni(gi,l,4)  =  < 


aiwiil  -  qi)  +  1, 

gi  <  1 

wi(l  -  gi)  +  1, 

1  <  gi  <  4 

(3iWi{4:  -  gi). 

4  <  gi 

q;2U'2(6  -  g2)  +  1, 

g2  <  6 

w^2(6  -  Qt)  +  1, 

6  <  g2  <  7 

/?2W2(7  -  g2), 

7  <  g2. 

«2(g2,6,  7)  =  < 


These  equations  are  used  by  Equation  (3.10),  the  single  objective  to  be  solved  by  the 
subproblem. 
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3. 2. 1.3  MGPS-RS  for  Problems  with  Linear  Constraints.  This  step 
uses  the  NOMADm  [3]  implementation  of  MGPS-RS  to  solve  each  single  objective 
subproblem  of  the  form  discussed  in  Section  3.1.4.  MGPS-RS  is  discussed  in  detail 
in  Section  3.1.2  and  has  been  shown  to  have  strong  convergence  properties.  (See 
Section  3.3.3. 1  and  [137].) 

3. 2. 1.4  MV-MADS  for  Problems  with  Nonlinear  Constraints.  Sim¬ 
ilarly,  this  step  of  SMOMADS  uses  the  NOMADm  implementation  of  MVMADS- 
RS  [3]  to  solve  each  single  objective  subproblem  of  the  form  discussed  in  Section 
3.1.4.  MADS  is  discussed  in  detail  in  Section  3.1.3  and  has  also  been  shown  to 
have  strong  convergence  properties  for  the  deterministic  case  [11].  (See  Section  3.3.1 
and  [5]  for  MV-MADS.)  Gonjecture  of  convergence  results  for  the  stochastic  case 
(MVMADS-RS)  are  discussed  in  Section  3.3.10. 

3.2. 1.5  Adding  Points  to  the  Efficient  Set.  The  algorithm  is  designed 
so  that  each  subproblem  should  produce  an  efficient  point.  In  fact,  this  is  always  the 
case  for  deterministic  problems  (see  Lemma  3.3.11).  In  stochastic  problems,  as  the 
number  of  iterations  of  the  single  objective  solver  is  allowed  to  approach  inhnity,  the 
solution  converges  to  an  efficient  point  with  probability  one  (see  Theorem  3.3.12  and 
Gonjecture  3.3.10).  However,  in  practice,  since  the  number  of  iterations  is  hnite,  the 
addition  of  dominated  points  is  possible.  Therefore,  points  are  checked  for  Pareto 
dominance  before  they  are  allowed  to  enter  the  efficient  set.  If  a  point  is  dominated, 
it  is  “hltered”  and  does  not  enter  the  efficient  set.  Ideally,  the  hlter  should  also 
check  to  see  if  the  new  point  dominates  other  points  in  the  current  efficient  set.  This 
feature  is  suggested  for  follow-on  research  in  Ghapter  6. 2. 3. 3. 

3. 2. 1.6  Tests  for  Quality  of  the  Pareto  Set.  An  exact  Pareto  set 
often  has  an  inhnite  number  of  efficient  points.  Since  multi-objective  solvers  provide 
only  an  approximate  Pareto  set,  the  quality  of  the  approximation  is  of  particular 
interest.  Relatively  few  papers  in  the  literature  focus  on  quality  metrics  for  Pareto 
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set  approximations  and  most  make  the  assumption  that  the  true  set  is  known  a  priori. 
Since  no  such  assumption  is  made  here,  the  quality  metrics  introduced  by  Wu  and 
Azarm  [163]  are  used  to  assess  the  quality  of  the  Pareto  set,  because  their  metrics 
measure  the  quality  (accuracy,  spread,  cluster,  etc.)  without  any  knowledge  of  the 
true  Pareto  set.  Specihc  metrics,  as  introduced  by  Wu  and  Azarm  [163],  are  described 
below: 

1.  Scaled  Objective  Space.  All  the  quality  metrics  used  in  this  research  are 
based  on  the  existence  of  “bad”  objective  function  values  {/(*,...,  /j}  that  are 
worse  than  all  points  in  the  Pareto  set,  and  “good”  objective  function  values 
{/f , . . . ,  /j}  that  are  better  than  all  points  in  the  Pareto  set.  The  utopia  and 
nadir  points  can  be  used  for  these  values,  other  points,  which  dehne  a  smaller 
region,  can  also  be  used,  as  determined  by  the  user.  All  the  metrics  require 
that  Pareto  optimal  points  be  scaled  using 

(3.13) 

■>  j 

where  fj{xk)  is  the  value  of  the  j-th  objective  at  the  point  Xk,  f  j{xk)  is  the 
scaled  value  of  the  j-th  objective  at  the  point  Xk,  and  /J  and  fj  are  the  good 
and  bad  values,  respectively,  of  the  j-th  objective  function.  The  metrics  that 
follow  are  heavily  dependent  on  choice  of  and  points.  Thus  these  metrics 
are  better  suited  to  relative  comparisons  between  sets,  e.g.,  measuring  how  the 
quality  of  the  set  improves  as  points  are  added,  rather  than  giving  an  absolute 
measure  of  the  quality  of  a  particular  set. 

2.  Hyperarea  Difference.  The  hyperarea  difference  {HD{P))  provides  a  quan¬ 
titative  means  to  evaluate  the  difference  between  the  size  of  the  scaled  objective 
space  dominated  by  an  approximate  Pareto  set  and  that  of  the  space  dominated 
by  the  true  set  itself  (which  is  equal  to  1  when  scaled).  Although  the  actual 
Pareto  set  is  usually  unknown,  it  is  still  possible  to  instead  compare  two  ex¬ 
perimental  sets  with  respect  to  how  much  worse  they  are  than  the  true  Pareto 
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in  the  scaled  objective  space.  The  salient  point  is  that,  thongh  the  sets  are 
measnred  with  respect  to  the  trne  set,  the  comparison  and  inferences  are  of 
the  two  experimental  sets  in  relation  to  each  other.  Given  an  observed  Pareto 
optimal  set,  scaled  appropriately  by  Eqnation  (3.13),  the  hyperarea  difference 
can  be  calcnlated  as 

HD{P)  =  1- 

rVp  m  V  _  1111 

2  n  1  -  max  (/i(a;fcj)  U 

kr=kr^H=^  L  ^  ^  J  J  J  J 

where  np  is  the  nnmber  of  observed  Pareto  optimal  points.  The  hyperarea 
metric  becomes  difficnlt  to  calcnlate  as  the  nnmber  of  points  in  the  efficient 
set  becomes  large.  For  this  reason,  this  metric  was  not  calcnlated.  However, 
it  conld  be  nsefnl  in  the  determination  of  termination  criteria,  as  described  in 
Section  4.4.2. 

3.  Pareto  Spread.  The  Pareto  spread  measnres  how  widely  the  experimental 
Pareto  optimal  set  spreads  over  the  objective  space.  The  overall  Pareto  spread 
{OS{P))  considers  all  objectives  together  and  is  defined  as  the  volnme  ratio  of 
two  hyper-rectangles.  One  rectangle  is  defined  by  and  for  each  objec¬ 
tive  and  the  other  is  defined  by  the  extreme  points  in  the  experimental  Pareto 
optimal  set.  The  overall  Pareto  spread  can  calcnlated  as 

m 

OS{P)  =  JJ  \max  {fi{xk)  :  k  =  l,2,.. .  ,np}  -  min  {fi{xk)  :  k  =  l,2,.. .  ,np}  \  . 

i=l 

Similarly,  the  k-th  objective  Pareto  spread 
OSk{P)  =  [max  [fk  (xi)  :i  =  l,2,...  ,np}  -  min  [fk  {xi)  :i  =  l,2,.. .  ,np}\ 
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is  a  measure  of  how  the  experimental  Pareto  optimal  set  spreads  over  the  ob¬ 
jective  space  when  each  objective  is  considered  individually. 

4.  Number  of  Distinct  Choices.  When  comparing  two  Pareto  sets,  the  one 
with  more  points  may  appear  to  be  better.  However,  this  may  not  be  the  case. 
If  the  points  in  a  Pareto  set  P  are  too  close  together,  according  to  some  distance 
metric,  they  may  then  be  indistinguishable  from  each  other  for  practical  design 
purposes.  The  number  of  distinct  choices  (NDCfj,{P))  is  a  measure  of  the 
number  of  Pareto  points  that  are  distinguishable  {i.e.,  far  enough  away)  from 
each  other.  It  is  calculated  as  follows.  A  quantity  yU  is  given  by  the  decision 
maker  which  divides  the  m-dimensional  objective  space  into  1/yu"*  small  grids. 
The  indicator  variable 


NT^{q,P) 


1,  3  Pfc  e  P  Pfc  e  T^{q) 
0,  y  Pk^  P  Pk^  T^{q) 


is  used  to  indicate  whether  a  particular  region  P^(g)  contains  any  point  pk  in  the 
Pareto  set.  The  number  of  distinct  choices  can  then  be  calculated  by  summing 
the  number  of  regions  that  contain  at  least  one  point, 

v—1  v—1  v—1 

l7n=0  l2=0ll=0 

where  q  =  (gi,  g2,  •  •  • ,  qm)  with  qi  =  -. 

V 

5.  Cluster.  Like  the  previous  metric,  the  cluster  metric  CL^[P)  quantihes  how 
closely  spaced  together  and,  as  can  be  seen  in 

CL  (P)  = 

^  '  NDC,,(P)’ 

it  is  the  ratio  of  the  number  of  experimental  Pareto  solutions  N(P)  to  the 
number  of  distinct  choices. 
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(a)  Objective  Space  (b)  Design  Space  Showing  Mesh 

Figure  3.8:  Notional  approximated  Pareto  set  and  mesh  for  a  two-objective  problem 
with  two  design  variables 


3.2.2  Stage  2.  Though  interactive  specihcation  of  aspiration/reservation 
levels  and  scalarization  functions  of  the  type  discussed  in  Section  3.1.4  is  capable  of 
hnding  Pareto  points  in  nonconvex  regions  of  the  Pareto  front,  it  is  not  guaranteed  to 
hnd  all  points.  It  is  possible  that  some  efficient  points  may  not  be  found.  Therefore, 
a  second  (optional)  stage  is  added  for  those  cases  in  which  missing  points  may  pose 
a  particular  problem.^  This  stage  does  not  effect  the  convergence  theory.  In  this 
stage,  a  discrete  mesh  (similar  to  that  used  for  MGPS-RS  and  MADS)  around  the 
current  efficient  points  is  generated  (see  Section  3.1.2)  and  then  a  multi-objective 
ranking  and  selection  algorithm  is  used  to  check  for  new  efficient  points  on  the  mesh, 
which  is  similar  to  the  frame  used  by  MADS  given  in  Equation  (3.4),  except  that  So  is 
replaced  by  Se,  the  set  of  efficient  points  found  in  stage  1.  A  graphical  representation 
of  a  notional  problem  is  shown  in  Figure  3.8. 


3.2.2. 1  Multi- Objective  Ranking  and  Selection.  A  version  of  Multi¬ 
objective  Optimal  Computing  Budget  Allocation  algorithm  (MOCBA),  developed  by 
Lee  et  al.  is  used  to  check  for  new  efficient  points  on  the  mesh.  As  discussed  in  Section 

^In  this  research,  the  algorithm  performed  so  well  on  the  test  problems  that  the  second  stage  was 
not  necessary;  however,  it  is  suggested  as  follow-on  research  in  Section  6.2.2. 
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3.1.5,  MOCBA  has  been  used  successfully  for  multi-objective  ranking  and  selection 
problems  [83, 86]  and  the  observed  Pareto  set  determined  by  MOCBA  approaches  the 
true  Pareto  set  asymptotically  with  probability  1  [85]. 

3. 3  Convergence  Results 

Since  existing  convergence  results  for  MADS  apply  only  to  NLP  problems,  the 
hrst  part  of  this  subsection  provides  a  convergence  theory  for  an  extension  of  MADS 
to  mixed  variable  problems,  called  MV-MADS.  The  extension  is  very  straightforward, 
so  that  the  algorithm  is  almost  identical  to  mixed  variable  pattern  search,  except  for 
the  traditional  differences  between  GPS  and  MADS  with  respect  to  poll  directions 
and  and  mesh  parameters.  This  work  is  based  on  the  joint  work  done  in  [5],  and  the 
results  presented  here  match  those  of  [5],  unless  otherwise  noted. 

3.3.1  Convergence  Results  for  Mixed-Variable  MADS  (MV-MADS).  Given 
a  constrained  optimization  problem  with  objective  function  /  :  D  — >  R,  where  D 
represents  the  feasible  region,  let  /q  :  R"  — >  R  be  dehned  by  /q  =  f  -\-  tf,  where 

=  0  if  a;  G  D,  and  -|-oo  otherwise,  and  assume  the  following: 

Al.  An  initial  point  Xq  with  ffi(xo)  <  oo  is  available. 

A2.  All  iterates  {xa,}  generated  by  MV-MADS  he  in  a  compact  set. 

A3.  The  set  of  discrete  neighbors  Af{xk)  lies  on  the  mesh  M^. 

Under  these  assumptions,  the  following  results  are  obtained  by  proofs  that  are  iden¬ 
tical  to  those  found  in  Audet  and  Dennis  [9]  and  Abramson  [2]  for  mixed  variable 
GPS: 


•  liminfA^  =  liminfA™  =  0; 

fc-H>-f-00  /c— ^H-OO 

•  there  exists  a  refining  subsequence  {xk}kGK  of  minimal  frame  centers  for  which 
there  are  limit  points  x  =  lim  Xk,  y  =  lim  yk,  and  5  =  (z'^,  y'^)  =  lim  Zk,  where 

k'^K  kGK  k^K 

each  G  G  is  the  endpoint  of  the  extended  poll  step  initiated  at  yk  G  jV{xk), 
and  lim  A?  =  0. 

keK 
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Some  of  the  results  that  follow  require  the  additional  assumption  that  y  G  N{x). 
Additionally,  note  the  following  properties  from  Clarke  [30]: 

•  Any  convex  set  is  regular  at  each  of  its  points. 

.  Both  Tg°(  x)  and  T^\x)  are  closed,  and  both  T^\x)  and  T^{x)  are  convex. 

.  Ti^{x)  C  TE\x)  C 

•  Rockafellar  [124]  showed  that,  if  T^{x)  is  nonempty,  T^{x)  =  mt{T^\x)) ,  and 
therefore,  T^\x)  =  d{T^{x)). 

A  generalization  [74]  of  the  Clarke  [30]  directional  derivative,  in  which  function  eval¬ 
uations  are  restricted  to  points  in  the  domain,  is  needed.  The  generalized  directional 
derivative  of  a  locally  Lipschitz  function  /  at  a;  G  hi  in  the  direction  v  G  is  dehned 
by 

rix^v)  :=  hmsup  “ /W .  (3,14) 

y  ^  X,  yen 
t  I  0,  y  +  tv  e  n 

The  following  four  dehnitions  [30, 74, 124]  are  needed.  They  have  been  adapted  to  the 
context  of  this  particular  problem,  where  only  a  subset  of  the  variables  are  continuous. 
The  standard  dehnitions  follow  when  all  variables  are  continuous,  he.,  when  n'^  =  0 
and  X  =  x^. 

Definition  3.3.1.  A  vector  v  G  M"’'"  is  said  to  he  a  hypertangent  vector  to  the  con¬ 
tinuous  variables  of  the  set  hi  at  the  point  x  =  {x^^,  x'^)  e  fl  if  there  exists  a  scalar 
e  >  0  such  that 

{y  tw,x^)  G  hi  for  all  y  G  B^{x^)  with  {y,x^)  e  fl,  w  e  B^{v)  and  0  <  t  <  e. 

The  set  Tq(x)  of  all  hypertangent  vectors  to  ft  at  x  is  called  the  hypertangent  cone 
to  hi  at  X. 
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Definition  3.3.2.  A  vector  v  G  M”’'"  is  said  to  be  a  Clarke  tangent  vector  to  the 
continuous  variables  of  the  set  at  the  point  x  =  (x‘^,  E  cl(li)  if  for  every  sequence 
{yk}  that  converges  to  x^  with  {yk,x'^)  E  and  for  every  sequence  of  positive  real 
numbers  {tk}  converging  to  zero,  there  exists  a  sequence  of  vectors  {wk}  converging 
to  V  such  that  {yk  +  tkWk,x‘^)  E  The  set  Tq\x)  of  all  Clarke  tangent  vectors  to  Q, 
at  X  is  called  the  Clarke  tangent  cone  to  ft  at  x. 

Definition  3.3.3.  A  vector  v  E  K"'"  is  said  to  be  a  tangent  vector  to  the  continuous 
variables  of  the  set  ft  at  the  point  x  =  {x‘^,x‘^)  E  cl(fi)  if  there  exists  a  sequence  {yk} 
that  converges  to  x^  with  (yk,x'^)  G  and  a  sequence  of  positive  real  numbers  {Afc} 
for  which  v  =  liuik  XkiVk  —  The  set  Tf{°{x)  of  all  tangent  vectors  to  ft  at  x  is 
called  the  contingent  cone  to  fl  at  x. 

Definition  3.3.4.  The  set  is  said  to  be  regnlar  at  x  ifTff^x)  =  TR\x). 

The  resnlts  of  this  section  make  nse  of  a  generalization  [74]  of  the  Clarke  [30] 
directional  derivative,  in  which  fnnction  evalnations  are  restricted  to  points  in  the 
domain.  Fnrthermore,  the  notions  of  generalized  directional  derivatives  and  gradient 
are  restricted  to  the  snbspace  of  continnons  variables.  The  generalized  directional 
derivative  of  a  locally  Lipschitz  fnnction  f  at  x  =  (x^,  x^)  G  hi  in  the  direction 
V  E  K”'"  is  dehned  by 

fof  N  r  f{y  +  tv,x 

f  [x-,v)  :=  hmsnp  - 

y  — >  (y,  x'*)  G  ST 
i  i  0,  (y  +  tv,xA  e  O 

Fnrthermore,  it  is  shown  in  [11]  that  if  Tff[x)  is  not  empty  and  v  E  T^\x),  then 

f°{x-,v)  =  lim  f°{x-,u).  (3.16) 

U  ^  V, 

ueT^ix) 

Other  notions  regarding  derivatives  are  generalized  as  follows.  Denote  by  V  f{x)  E 
R"'"  and  df{x)  C  R”'",  respectively,  the  gradient  and  generalized  gradient  of  the  fnnc- 


)  - 


(3.15) 
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tion  f  at  X  =  G  Vt  with  respect  to  the  continuous  variables  x^  while  holding 

the  categorical  variables  x'^  constant.  In  particular,  the  generalized  gradient  of  /  at 
X*  [30]  with  respect  to  the  continuous  variables  is  dehned  by 

df{x)  :=  {s  e  M"'''  :  f°{x]v)  >  v^s  for  all  v  G  . 

The  function  /  is  said  to  be  strictly  differentiable  at  x  with  respect  to  the 
continuous  variables  if  the  generalized  gradient  of  /  with  respect  to  the  continuous 
variables  at  a;  is  a  singleton;  he.,  df(x)  =  {V/(a;)}. 

The  hnal  dehnition,  which  is  adapted  from  the  original  MADS  algorithm  in 
[11]  for  the  mixed  variable  case,  provides  some  terminology  for  stationarity  in  the 
nonsmooth  case. 

Definition  3.3.5.  Let  f  be  Lipschitz  near  x*  G  D.  Then  x*  is  said  to  be  a  Clarke, 
or  contingent,  stationary  point  of  f  over  ft  with  respect  to  the  continuous  variables  if 
f°{x*;  v)  >0  for  every  direction  v  in  the  Clarke  tangent  cone  to  D  at  x* ,  or  contingent 
cone  to  D  at  x* ,  respectively. 

In  addition,  x*  is  said  to  be  a  Clarke,  or  contingent,  KKT  stationary  point  of 
f  over  D  if  —Vf(x*)  exists  and  belongs  to  the  polar  of  the  Clarke  tangent  cone  to  kl 
at  X* ,  or  contingent  cone  to  ft  at  x* ,  respectively. 

If  Ll'^{x‘^)  =  R"'"  or  x*^^  lies  in  the  relative  interior  of  Ct^{x'^),  then  a  stationary 
point  as  described  by  Dehnition  3.3.5  meets  the  condition  that  f°{x*]v)  >  0  for  all 
V  G  R”'".  This  is  equivalent  to  0  G  df{x*). 

Main  convergence  results  for  the  MVMADS  algorithm  consist  of  four  theorems, 
all  of  which  are  generalizations  of  similar  results  from  MADS  [11]  or  mixed  variable 
pattern  search  [2,6,9].  The  hrst  result  establishes  a  notion  of  directional  stationarity 
at  certain  limit  points,  and  the  second  ensures  local  optimality  with  respect  to  the  set 
of  discrete  neighbors.  The  remaining  two  results  establish  Clarke-based  stationarity 
in  a  mixed  variable  sense. 
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Theorem  3.3.6.  Let  w  be  the  limit  point  of  a  refining  subsequence  or  the  associ¬ 
ated  subsequence  of  extended  poll  endpoints,  and  let  v  be  a  refining  direction  in 
the  hypertangent  cone  Tff  (w)  ■  If  f  is  Lipschitz  at  w  with  respect  to  the  continuous 
variables,  then  f°(w;  (t,0))  >  0. 


Proof.  Let  K  be  an  indexing  set  of  {wk\k(^K-,  a  refining  snbseqnence  converging  to 
w  =  {w^,w^).  Withont  any  loss  of  generality,  assnme  that  Wk  =  (wl,w^)  for  all 
k  ^  K.  In  accordance  with  Dehnition  3.2  in  [11],  let  v  =  liuik^Lj^^  £  'I'ni'^)  be  a 
rehning  direction  for  w,  where  dk  G  Dk  for  all  k  ^  L  and  L  is  some  snbset  of  K. 

Becanse  Wk  converges  to  w  and  by  the  dehnition  of  the  MADS  poll  set,  A™||(ifc|| 
is  bonnded  above  by  max{||d'||  :  d'  G  D}  where  D  is  a  hnite  set  of  directions  and 
Ap  converges  to  zero,  it  follows  that  A™||dfc||  also  converges  to  zero.  Thns,  it  also 
follows  from  Eqnation  (3.15)  that: 


f°{w;v)  =  limsnp 

y  — »  w'^,  {y,  w'^)  €  O 
MO,  (y  +  tu,  w’^)  G  Q 
u  ^  v,u  G  (w) 


f{y  +  tu,w^)  -  f{y,w'^) 


>  lim  snp 
k  e  L 


/  (<  + A™||4||^,w‘^)  -f{wl,w^) 


Arii4ii 


f  {wk  +  Al‘dt,‘w'‘)  -  f{wt) 

=  limsnp  — ^ - r-rrrwr-rv -  >  U. 


k  e  L 


The  last  ineqnality  holds  becanse  (wk  +  A^dk,  vL^)  G  D  and  f{wk  +  A^dk,  vL^)  >  f{wk) 
(since  Wk  is  a  minimal  frame  center)  for  all  snfficiently  large  k  &  L.  □ 

The  next  resnlt  gives  snfficient  conditions  nnder  which  a;  is  a  local  minimizer 
with  respect  to  its  discrete  neighbors. 

Theorem  3.3.7.  If  f  is  lower  semi- continuous  at  x  and  upper  semi-continuous  at 
y  G  M{x)  with  respect  to  the  continuous  variables,  then  f(x)  <  fifi). 
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Proof.  Since  k  &  K  ensures  that  {xk}k&K  are  minimal  frame  centers,  it  follows  that 
f{xk)  <  fivk)  for  all  k  ^  K.  By  the  assumptions  of  lower  and  upper  semi-continuity 
on  /  and  the  dehnitions  of  x  and  it  also  follows  that  f(x)  <  lim^gx /(a^fc)  < 

\imk(.K  fiVk)  =  fiy)-  □ 

The  next  theorem  lists  conditions  that  ensure  that  x  satishes  certain  stationary 
conditions,  under  various  smoothness  requirements. 

Theorem  3.3.8.  Assume  that  Tff{x)  ^  0  and  the  set  of  refining  directions  is  asymp¬ 
totically  dense  in  Tff{x). 

1.  If  f  is  Lipschitz  near  x  with  respect  to  the  continuous  variables,  then  x  is  a 
Clarke  stationary  point  of  f  on  ft  with  respect  to  the  continuous  variables. 

2.  If  f  is  strictly  differentiable  at  x  with  respect  to  the  continuous  variables,  then 
X  is  a  Clarke  KKT  stationary  point  of  f  on  ft  with  respect  to  the  continuous 
variables. 

Furthermore,  if  hi  is  regular  at  x,  then  the  following  hold: 

1.  If  f  is  Lipschitz  near  x  with  respect  to  the  continuous  variables,  then  x  is  a 
contingent  stationary  point  of  f  on  ft  with  respect  to  the  continuous  variables. 

2.  If  f  is  strictly  differentiable  at  x  with  respect  to  the  continuous  variables,  then 
X  is  a  contingent  KKT  stationary  point  of  f  on  fl. 

Proof.  First,  Rockafellar  [124]  showed  that  if  the  hypertangent  cone  is  not  empty  at 
X,  then  Tq\x)  =  cl{Tff{x)).  Since  the  set  S  of  rehning  directions  for  /  at  a;  is  a  dense 
subset  of  Tq{x),  S  is  also  a  dense  subset  of  Tf^\x).  Thus,  any  vector  v  G  Tf{\x)  can 
be  expressed  as  the  limit  of  directions  in  S,  and  the  hrst  result  follows  directly  from 
(3.16)  and  Theorem  3.3.6. 

Strict  differentiability  ensures  the  existence  of  Vf{x)  and  that  Vf{x)'^v  = 
f°{x;v)  for  all  v  G  T^\x).  Since  f°{x;v)  >  0  for  all  v  G  T^fix),  it  follows  that 
(— V/(a;))^n  <  0,  and  the  second  result  follows  from  Definition  3.3.5.  Furthermore, 
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if  is  regular  at  x,  then  by  Definition  3.3.4,  T^\x)  =  Tq°{x),  and  the  hnal  two 
results  follow  directly  from  Dehnition  3.3.5.  □ 

The  next  result  is  similar  to  Theorem  3.3.8  but  considers  the  limit  of  extended 
POLL  endpoints  5  instead  of  x. 

Theorem  3.3.9.  Assume  thaty  G  M(x),  T^{z)  ^  0,  and  the  set  of  refining  directions 
is  asymptotically  dense  in  Tff{z). 

1.  If  f  is  Lipschitz  near  z  with  respect  to  the  continuous  variables,  then  z  is  a 
Clarke  stationary  point  of  f  on  ft  with  respect  to  the  continuous  variables. 

2.  If  f  is  strictly  differentiable  at  z  with  respect  to  the  continuous  variables,  then 
z  is  a  Clarke  KKT  stationary  point  of  f  on  ft  with  respect  to  the  continuous 
variables. 

Furthermore,  if  D  is  regular  at  z,  then  the  following  hold: 

1.  If  f  is  Lipschitz  near  z  with  respect  to  the  continuous  variables,  then  z  is  a 
contingent  stationary  point  of  f  on  ft  with  respect  to  the  continuous  variables. 

2.  If  f  is  strictly  differentiable  at  z  with  respect  to  the  continuous  variables,  then 
z  is  a  contingent  KKT  stationary  point  of  f  on  D. 

Proof.  The  proof  is  identical  to  that  of  Theorem  3.3.8,  but  with  5  replacing  x.  □ 

3.3.2  Convergence  of  MV-MADS  with  Ranking  and  Selection  (MVMADS-RS). 
Convergence  for  mixed  variable  MADS  for  the  stochastic  case  has  not  yet  been 
formally  proven.  Because  of  the  similarities  between  MV-MADS  and  MGPS-RS 
discussed  in  Section  3.1.3  and  the  rigorous  convergence  results  that  exist  for  MGPS- 
RS,  the  following  conjecture  is  made  with  the  belief  that  a  proof  similar  to  that  for 
MGPS-RS  would  follow  for  MVMADS-RS. 

Conjecture  3.3.10.  Suppose  the  sequence  of  iterates  generated  by  MVMADS-RS 
converges  to  x  E  D.  Then  x  meets  the  first-order  necessary  conditions  (in  the  forms 
listed  below)  for  optimality  a.s.: 
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•  If  f  is  Lipschitz  near  x,  then  x  is  a  Clarke  stationary  point  of  f  onVL. 

•  If  f  is  strictly  differentiable  at  x  and  Tff[x)  7^  0,  then  x  is  a  Clarke  KKT 
stationary  point  of  f  over  ft . 

•  If  f  is  strictly  differentiable  at  x,  ft  is  regular  at  x,  and  Tff[x)  7^  0,  then  x  is  a 
contingent  KKT  stationary  point  of  f  over 

3. 3. 3  Convergence  Results  for  Multi- Objective  Problem.  Convergence  resnlts 
for  the  mnlti-objective  problem  are  now  presented.  First,  it  is  shown  that  if  the 
snbproblems  converge  to  an  optimal  solntion,  that  solntion  is  also  Pareto  optimal. 

Lemma  3.3.11.  Civen  a  feasible  point  x*  E  VL  of  the  stochastic,  multi-objective, 
mixed-variable  optimization  problem  defined  in  Equation  (1.5),  if  x*  is  a  global  mini- 
mizer  of  a  convex  combination  of  the  J  objectives,  then  x*  is  Pareto  optimal. 

Proof.  Assnme  to  the  contrary  that  x*  is  not  Pareto  optimal.  If  x*  is  not  Pareto 
optimal,  by  Definition  1.1.3,  there  exists  some  x  E  D  snch  that  Fk(x)  <  Fk(x*) 
for  /c  =  1,2,...,J  and  Fi{x)  <  Fi{x*)  for  some  i  G  {1,2,...,J}.  Thns,  the 

j  j 

positive  snm,  '^CiFfix)  <  '^ciFfix*),  which  contradicts  the  assnmption  that  x*  = 

i=l  i=l 

argmin  j  '^CiFi(x)  j.  Therefore,  x*  is  Pareto  optimal.  □ 

xGfi  V*=l  / 

3.3.3. 1  Linearly  Constrained  Problems.  Convergence  resnlts  for  SMO- 
MADS  depend  entirely  on  the  resnlts  for  the  snb-problems.  The  following  assnmptions 
are  made  [139]: 

1.  All  iterates  prodnced  by  the  the  MGPS-RS  algorithm  he  in  a  compact  set. 

2.  The  objective  fnnction  /  is  continnonsly  differentiable  with  respect  to  the  con- 
tinnons  variables  when  the  discrete  variables  are  fixed. 

3.  For  each  set  of  discrete  variables  X‘^,  the  corresponding  set  of  directions  H*  = 
GiZi  inclndes  tangent  cone  generators  for  every  point  in 

4.  The  rnle  for  selecting  directions  D],  conforms  to  for  some  £  >  0. 
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5.  For  each  q  =  1,2, .. .  ,nc,  the  responses  are  independent,  identically 

and  normally  distribnted  random  variables  with  mean  /  (X^)  and  nnknown 
variance  cr^  <  oo,  where  af  ^  whenever  I  ^  q. 

6.  For  the  r-th  R&S  procednre  with  candidate  set  C  =  {Yi,  Y2, . . . ,  Yn^},  RS  {C,  <5^) 

gnarantees  correctly  selecting  the  best  candidate  Y^i]  G  C  with  probability  of  at 
least  1  —  whenever  /  (F[g])  —  /  (F[i])  >  5r  for  any  g  G  {2,3,...,  nc}- 

7.  For  all  bnt  a  hnite  nnmber  of  MGPS-RS  iterations  and  snb-iterations,  the  best 
solntion  Yfij  G  (7  is  nniqne,  he.,  f  {Y[l])  7^  f  {^[q])  for  fol  Q  ^  (2,  3, . . . ,  nc} 
where  C  =  (Fi,  Y2, C  Mk  at  iteration  k. 

Given  a  feasible  point  a;*  G  G  of  a  stochastic,  mnlti-objective,  mixed-variable  opti¬ 
mization  problem  as  dehned  in  Eqnation  (1.5),  where  the  continnons  variables  in  Vt 
are  restricted  by  bonnd  and  linear  constraints,  the  following  key  resnlt  applies. 

Theorem  3.3.12.  The  sequence  of  iterates  generated  by  each  subproblem  of  stochastic 
multi-objective  pattern  search  (SMOMADS)  contains  a  limit  point  that  meets  the  first- 
order  necessary  conditions  for  Pareto  optimality,  almost  surely  (a.s.). 


Proof.  The  SMOMADS  algorithm  generates  each  snbproblem  as  a  nonnegative  com¬ 
bination  of  the  J  objectives  of  the  original  problem,  he.,  Z{x)  =  (^^CiFi(a;)Y  a  >  0. 


\*=i  / 

Each  snbproblem  is  then  solved  nsing  MGPS-RS.  Resnlts  from  Theorem  4.8  and  Theo¬ 
rem  4.13  [139]  show  that  the  seqnence  of  iterates  prodnced  in  the  snbproblem  contains 
a  limit  point  x*  satisfying  hrst-order  conditions  for  optimality  a.s.  By  Lemma  3.3.11, 
if  X*  is  optimal,  it  is  also  Pareto  optimal.  □ 


3. 3. 3. 2  Nonlinearly  Constrained  Problems.  The  following  assnmptions 
are  made  [5]: 

1.  All  iterates  Xk  prodnced  by  the  the  MVMADS-RS  algorithm  he  in  a  compact 
set. 

2.  The  objective  fnnction  /  is  continnonsly  differentiable  with  respect  to  the  con¬ 
tinnons  variables  when  the  discrete  variables  are  fixed. 
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3.  For  each  q  =  1,2, .. .  ,nc,  the  responses  are  independent,  identically 

and  normally  distribnted  random  variables  with  mean  /  (X^)  and  nnknown 
variance  cr^  <  oo,  where  af  ^  whenever  I  ^  q. 

4.  For  the  r-th  R&S  procednre  with  candidate  set  C  =  {Yi,  Y2, . . . ,  Yn^},  RS  {C,  ar,  Sr) 
gnarantees  correctly  selecting  the  best  candidate  Y^i]  G  C  with  probability  of  at 
least  1  —  ar  whenever  /  (F[g])  —  /  (h[i])  >  Sr  for  any  g  G  {2,3,...,  nc}- 

5.  For  all  bnt  a  hnite  nnmber  of  MVMADS-RS  iterations  and  snb-iterations,  the 

best  solntion  Yji]  G  C*  is  nniqne,  he.,  /  (R[l])  7^  /  (^<7])  for  all  g  G  (2,  3, ... ,  nc} 
where  C  =  (Yi,  12,  •  •  • ,  ^  at  iteration  k. 

Theorem  3.3.13.  Suppose  the  sequence  of  iterates  generated  by  a  subproblem  of 
SMOMADS  solved  using  MVMADS-RS  converges  to  x  E  Q.  Then  x  meets  the  first- 
order  necessary  conditions  (in  the  forms  listed)  for  optimality  a.s.: 

•  If  f  is  Lipschitz  near  x,  then  x  is  a  Clarke  stationary  point  of  f  on  D 

•  If  f  is  strictly  differentiable  at  x  and  Tff{x)  7^  0,  then  x  is  a  Clarke  KKT 
stationary  point  of  f  over  . 

•  If  f  is  strictly  differentiable  at  x,  D  is  regular  at  x,  and  Tff[x)  7^  0,  then  x  is  a 
contingent  KKT  stationary  point  of  f  over  fl. 

Further,  if  x  is  in  fact  optimal,  it  is  also  Pareto  optimal. 


Proof.  The  SMOMADS  algorithm  generates  each  snbproblem  as  a  nonnegative  com¬ 
bination  of  the  J  objectives  of  the  original  problem,  he.,  Z{x)  =  a  >  0. 


V*=i  / 

Each  snbproblem  is  then  solved  nsing  MADS.  Thns,  by  Conjectnre  3.3.10,  the  limit 
point  X  satishes  hrst-order  necessary  conditions  for  optimality,  he.,  is  a  stationary 
point,  a.s.  Therefore,  by  Lemma  3.3.11,  if  x  is  optimal,  it  is  also  Pareto  optimal.  □ 


3.4  Summary 

In  this  chapter,  the  elements  of  the  SMOMADS  algorithm — R&S,  MGPS-RS, 
MADS,  and  aspiration/reservation  level  analysis — have  been  described  and  the  algo- 
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rithm  itself  has  been  outlined  and  described.  Additionally,  it  has  been  shown  that 
each  subproblem  of  the  SMOMADS  algorithm  converges  almost  surely  to  stationary 
points  appropriately  dehned  in  the  mixed  variable  domain,  and  that  if  such  points 
are  globally  optimal,  then  they  are  also  Pareto  optimal.  After  the  methodology 
was  developed  and  theoretical  convergence  properties  determined,  algorithmic  imple¬ 
mentations  were  developed  for  testing.  These  implementations  are  given  in  Chapter 
IV. 
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IV.  Algorithmic  Implementation 

In  this  chapter,  the  details  of  the  various  SMOMADS  algorithm  implementations  are 
presented.  The  implementations  can  be  categorized  based  on  four  central  issues: 
determining  aspiration/reservation  level  pairings  in  the  objective  space  (he.,  in  the 
master  problem),  type  of  stochastic  solver,  type  of  solution  filter,  and  termination 
criteria.  These  issues  are  discussed  in  Sections  4. 1-4.4,  respectively.  The  implemen¬ 
tations  are  described  in  Section  4.5. 

4-1  Objective  Space  Search  Pattern 

As  discussed  in  Section  3. 2. 1.1,  specihc  choices  of  aspiration  and  reservation 
levels  are  used  to  combine  the  multiple  objectives  of  the  original  problem  into  a  single 
objective  subproblem,  which  generates  a  point  on  the  approximate  Pareto  front  in  a 
given  region  of  interest.  In  generating  an  approximate  Pareto  front,  it  is  important 
that  efficient  points  not  be  too  clustered  in  the  same  region,  so  that  the  true  Pareto 
front  can  be  more  fully  and  accurately  characterized.  Principles  of  experimental 
design  were  used  in  the  objective  space  to  accomplish  this  goal.  Specihcally,  three 
methods  were  chosen,  and  in  each  case,  the  levels  of  the  design  were  chosen  to  be 
specihc  predetermined  values  for  the  aspiration  and  reservation  levels. 

1.  Full  Factorial  Design.  The  full  factorial  design  (FFD)  has  as  design  points 
every  possible  combination  of  preselected  design  variable  values.  Though  full 
factorial  designs  provide  information  about  linear,  interaction,  and  quadratic 
effects  {i.e.,  in  characterizing  the  shape  of  the  Pareto  front),  they  become  im- 
practically  large  for  relatively  few  numbers  of  design  variables  and  levels  because 
the  number  of  design  variables  grows  twice  as  fast  as  the  number  of  objective 
functions. 

2.  Central  Composite  Design.  The  central  composite  design  (CCD)  is  a  hve- 
level  variance-optimal  design  used  to  ht  second  order  models.  It  is  considered 
by  proponents  of  experimental  design  [104]  to  be  quite  useful  for  sequential 
experimentation,  in  which  lower  cost  screening  experiments  are  hrst  performed 
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to  determine  regions  where  further  scrutiny  is  required  [104],  Information  about 
shape  of  the  Pareto  front  can  be  determined  with  relatively  few  design  points 
[104], 

3.  Box-Behnken  Design.  The  Box-Behnken  design  (BBD)  was  developed  as  a 
three-level  alternative  to  the  CCD.  It  is  a  spherical  design  that  provides  good 
coverage  of  the  design  space  in  general.  However,  because  it  is  spherical,  instead 
of  cuboidal,  it  should  not  be  used  if  the  decision  maker  is  particularly  concerned 
with  the  extreme  points  of  the  given  domain  [25, 104]. 

Though  only  these  three  methods  were  implemented  in  this  research,  other  experi¬ 
mental  designs  are  available  and  may  be  of  use.  Additional  insights  about  the  master 
problem  (he.,  the  Pareto  front)  may  allow  the  user  to  more  parsimoniously  allocate 
potentially  costly  objective  function  evaluations  (or  simulations).  Further  research 
in  this  area  is  suggested  in  Section  6.2.3. 

J^.2  Stochastic  Solvers 

Once  the  objectives  of  the  master  problem  are  combined  into  a  single  objective, 
the  MGPS-RS  and  MVMADS-RS  stochastic  optimization  algorithms,  discussed  in 
Sections  3.1.2  and  3.2. 1.4,  respectively,  were  applied  because  of  the  their  attractive 
convergence  properties  and  because  they  are  useful  on  industrial  problems  in  which 
derivatives  are  not  available  {e.g.  see  [144],  [17]  or  [4]).  The  actual  implementation 
makes  use  of  the  NOMADm  [3]  MATLAB®  software,  which  has  been  applied  suc¬ 
cessfully  to  several  single-objective  stochastic  test  problems  [45]  and  a  multi-echelon 
repair  system  optimization  application  [144]. 

NOMADm  was  extended  to  include  the  multi-objective  case  in  the  following 
way.  Code  for  interactive  specihcation  of  aspiration/reservation  level  analysis  was 
required  to  control  the  master  problem.  Since  no  suitable  MATLAB  code  exists 
for  doing  this,  new  code,  called  MOMADS  (see  Appendix  A),  was  written,  which 
incorporates  the  NOMADm  software  as  the  single-objective  stochastic  subproblem 
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solver  as  described  in  Section  3.2.  A  copy  of  this  new  code  is  inclnded  as  Appendix 

A. 

4-3  Solution  Filter 

For  deterministic  optimization  problems,  the  solntions  to  the  snbproblems  of 
the  SMOMADS  algorithm  are  gnaranteed  to  be  efficient,  or  non-dominated,  points 
(see  Section  3.3.3).  This  is  also  trne  of  stochastic  optimization  problems,  bnt  only  in 
the  limit,  he.,  if  the  snbproblem  is  allowed  an  inhnite  nnmber  of  iterations.  However, 
since  an  inhnite  nnmber  of  iterations  is  not  possible  in  practice,  some  non-efficient 
solntions  may  enter  the  efficient  set.  For  this  reason,  each  point  is  checked  for  Pareto 
dominance  before  entering  the  efficient  set,  and  if  it  is  dominated,  it  is  “hltered”  and 
does  not  enter  the  efficient  set. 

4-4  Algorithm  Termination  Criteria 

In  SMOMADS,  algorithmic  termination  mnst  occnr  both  in  the  snbproblem 
and  in  the  master  problem.  Termination  criteria  for  the  snbproblems  is  discnssed  in 
Section  4.4.1  and  termination  criteria  for  the  master  problem  is  discnssed  in  Section 

4.4.2. 


4-4  A  Terminating  the  Subproblems.  Even  if  an  algorithm  is  known  to  con¬ 
verge,  the  reality  of  imprecision  and  ronndoff  error  make  it  necessary  to  predeter¬ 
mine  stopping  criteria.  In  pattern  search  methods,  this  is  accomplished  by  stop¬ 
ping  the  algorithm  when  the  mesh  size  A™  falls  below  a  threshold  valne  A^;  i.e., 
Afc  <  At  [67,137].  This  criterion  has  been  nsed  extensively  in  practice,  and  Dolan  et 
al.  [44]  showed  that  for  problems  withont  noise  A™  provides  a  measnre  of  hrst-order 
stationarity;  he.,  ||V/  [x^)  jj  <  CA™,  where  C*  >  0  is  an  nnknown  constant.  Similar 
termination  criteria  may  be  nsed  for  MADS. 

In  a  stochastic  environment,  termination  criteria  are  typically  more  complex. 
Too  small  a  valne  of  At  may  increase  the  reqnired  nnmber  of  fnnction  evalnations 
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required  by  the  ranking  and  selection  portion  of  the  algorithm  to  an  unacceptable 
level,  whereas,  too  large  a  value  may  induce  premature  termination.  For  MGPS-RS 
applied  to  single  objective  problems  two  additional  termination  criteria  are  employed 
[137].  The  hrst  is  that  the  ratio  between  the  response  standard  deviation  S  and  the 
indifference  zone  parameter  Sr  exceeds  some  threshold;  namely, 

f  >  (4.1) 

where  K  may  be  selected  by  the  user,  based  on  available  sampling  budget;  larger 
values  of  K  allow  larger  budgeting  thresholds  [137].  This  criterion  is  a  measure  of 
signal-to-noise.  It  is  a  heuristic  for  signalling  that  a  high  response  noise  (compared 
to  the  indifference  zone  parameter)  is  causing  excessively  large  sampling  requirements 
[137] .  The  second  criterion  is  that  the  significance  level  in  the  ranking  and  selection 
procedure  satisfy 

(Xr  ^  Olj'-i 

where  ax  is  a  threshold  of  error  to  ensure  that  the  probability  of  correct  selection 
1  —  ut  from  among  the  candidates  is  sufficiently  high.  This  is  intended  to  prevent 
a  sequence  of  erroneous  selections  from  causing  the  step  size  to  decrease  to  a  small 
value  thereby  prematurely  terminating  the  algorithm. 

4-4-2  Terminating  the  Master  Problem.  In  Section  4.1  experimental  designs 
were  developed  to  search  the  objective  space.  In  those  implementations,  a  fixed 
number  of  design  points  are  determined  a  priori  and  the  algorithm  terminates  when 
these  experimental  design  points  are  exhausted;  however,  in  some  cases  this  may  not 
be  the  ideal  method  for  terminating  the  master  problem.  For  example,  in  the  full 
factorial  design,  the  number  of  design  points  can  be  prohibitively  large.  If  a  metric 
to  determine  the  quality  of  the  Pareto  set  existed,  it  could  be  used  to  terminate  the 
algorithm  without  running  all  the  design  points  if  the  Pareto  set  was  determined  to 
be  “adequate”.  Therefore,  one  area  of  future  work  is  to  apply  the  quality  metrics 
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introduced  by  Wu  and  Azarm  [163]  to  assess  the  Pareto  set. 


Figure  4.1:  Master  problem  algorithm  with  termination  check  (graphical) 


Termination  in  the  master  problem  was  accomplished  via  the  algorithm  depicted 
in  Figure  4.1  and  described  in  Figure  4.2.  The  approach  is  not  ideal,  as  it  requires 
intuitive  observation  on  the  part  of  an  experienced  analyst  with  moderate  insight 
into  the  expected  objective  space;  however,  for  the  purpose  of  this  initial  study,  it 
does  perform  quite  well.  Additionally,  in  many  real-world  engineering  problems, 
extreme  cases  are  not  typically  considered;  e.g.,  an  aircraft  designed  to  have  zero 
reliability  or  excessively  high  cost.  Furthermore,  some  information  about  the  decision 
maker’s  region  of  interest  is  likely  to  exist  and  provide  insight  into  whether  to  expend 
(objective  function  evaluation)  resources  to  hll  gaps  in  regions  of  little  or  no  interest 
to  the  decision  maker.  So,  though  not  theoretically  ideal,  this  approach  is  usually 
adequate  in  practice. 

It  has  been  suggested  previously  that  the  quality  metrics,  applied  a  posteriori 
in  this  study,  could  also  be  applied  after  each  successive  design  point.  This  would 
provide  insight  into  the  quality  of  the  Pareto  set  approximation  and  quantitative 
termination  criteria.  However,  this  would  only  provide  information  about  when  to 
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A  Method  For  Determining  When  to  Terminate  the  Master  Problem 

1.  Choose  initial  SMOMADS  implementation  based  on  problem  structure 
(types  of  constraints,  types  and  number  of  objectives,  etc.). 

2.  Run  the  implementation  chosen  in  step  1  and  graph  the  results  in  the  objec¬ 
tive  space.  (For  problems  with  more  than  three  objectives,  the  objectives 
can  be  graphed  pairwise.) 

3.  Visually  inspect  the  approximated  Pareto  front  and  consider  the  following: 

•  Do  gaps  appear? 

•  Is  the  spread  adequate  for  the  desired  range  of  aspiration  and  reserva¬ 
tion  levels? 

•  Are  desired  “extreme  solutions  ”  represented,  he.,  those  solutions  close 
in  one  objective  to  its  component  in  the  utopia  point? 

•  Are  solutions  clustered? 

4.  If  approximation  is  adequate,  then  keep  current  efficient  set  and  run  quality 
metrics.  Retain  current  Pareto  set  and  stop. 

5.  Otherwise,  determine  region  requiring  further  scrutiny,  e.g.,  regions  between 
points  A1  and  R1  and  between  points  A2  and  R2  in  Figure  4.3. 

Figure  4.2:  Master  problem  algorithm  with  termination  check 
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Figure  4.3:  Examining  missing  regions  of  the  Pareto  front 

terminate.  If  the  set  is  determined  to  be  inadequate,  the  metrics  by  themselves 
would  not  indicate  how  to  set  the  aspiration  and  reservation  levels  to  improve  the 
approximation  in  successive  runs.  The  analyst  would  still  be  needed  to  determine 
appropriate  settings. 

4-5  SMOMADS  Implementations 

Seven  instances  of  the  SMOMADS  algorithm,  dehned  by  the  choices  discussed  in 
Sections  4. 1-4.4  were  devised  and  tested.  They  are  listed  in  Table  4.1,  where  columns 
2-4  describe  the  choice  that  dehne  the  number  listed  in  column  1. 

4-6  Testing  the  Algorithm 

In  this  chapter,  implementation  of  the  SMOMADS  algorithm  has  been  described 
in  detail.  In  particular,  four  central  issues  were  addressed:  search  pattern  in  the 


72 


Nnmber 

Table  4.1: 
Solver 

SMOMADS  Implementations  Tested 
Experimental  Design  Filter  Test  Problems 

Section 

1 

MGPS-RS 

Fnll  Factorial 

No 

Viennet  4 

5. 1.1.1 

2 

MGPS-RS 

Fnll  Factorial 

Yes 

Viennet  3 

5. 1.1.2 

Poloni 

5. 1.1. 3 

Dias  PI 

5. 1.2.1 

Dias  r2 

5. 1.2.2 

Fonseca  FI 

5. 1.2.3 

Schaffer  F3 

5. 1.2.4 

DTLZ7 

5. 1.2.6 

3 

MGPS-RS 

Central  Composite 

No 

Viennet  4 

5. 1.1.1 

4 

MGPS-RS 

Central  Composite 

Yes 

Dias  PI 

5. 1.2.1 

Fonseca  FI 

5. 1.2.3 

Schaffer  F3 

5. 1.2.4 

DTLZ7 

5. 1.2.6 

5 

MGPS-RS 

Box-Behnken 

No 

Viennet  4 

5. 1.1.1 

6 

MVMADS-RS 

Fnll  Factorial 

Yes 

Tamaki 

5. 1.1.4 

Srinivas 

5. 1.2.5 

7 

MVMADS-RS 

Central  Composite 

Yes 

Tamaki 

5. 1.1.4 

Srinivas 

5. 1.2.5 

Disk  Brake 

5. 1.3.1 

objective  space  (be.,  in  the  master  problem),  type  of  stochastic  solver,  termination 
criteria,  and  whether  or  not  a  solntion  hlter  (as  described  in  Section  3. 2. 1.5)  was  nsed. 
The  seven  implementations  shown  in  Table  4.1  were  nsed  to  solve  pnblished  mnlti- 
objective  test  problems  with  known  resnlts  and  then  applied  to  a  mnlti-objective 
design  optimization  problem.  The  test  problems,  application  problem,  and  resnlts 
are  given  in  Chapter  V. 
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V.  Computation  and  Evaluation 

The  accuracy  of  the  SMOMADS  algorithm  was  verihed  by  its  use  on  a  set  of  test  prob¬ 
lems  with  known  solutions.  These  test  problems  were  modihed  (via  the  MATLAB® 
rand  function)  to  introduce  random  noise  into  each  objective  function  to  simulate  the 
algorithm’s  use  on  a  stochastic  system.  The  problems  (except  the  airfoil  design  prob¬ 
lem)  were  solved  in  Matlab  7.2.0  on  a  2.13  GHz  Pentium(R)M  processor  with  1GB  of 
RAM.  Initially,  a  prototype  of  the  MOMADS  code  (see  Appendix  A)  was  verihed  by 
its  use  on  one  test  problem  [149]  before  checking  an  entire  set  of  problems.  Results 
of  the  initial  test,  which  are  also  included  in  [150],  are  shown  in  Section  5. 1.1.1. 

Following  positive  initial  results,  additional  test  problems  were  identihed  target¬ 
ing  several  characteristics:  bi-objective,  multi-objective,  continuous-variable,  mixed- 
variable,  discontinuous  Pareto  front.  Specihc  test  problems  are  listed  in  Table  5.1. 
Results  for  each  test  problem  are  shown  in  Section  5.1  and  associated  quality  metrics 
are  listed  in  Table  5.3  in  Section  5.3. 


Table  5.1:  Published  Problems  Tested  using  SMOMADS 


Test  Problem 

#  Var 

#Obj 

Type  of 
Variables 

Sequential 

Experiment 

^  Test 
Points 

Results 
in  Section 

Viennet  4 

2 

3 

Gontinuous 

No 

4,209 

5.1. 1.1 

Viennet  3 

2 

3 

Gontinuous 

No 

4,096 

5. 1.1.2 

Poloni 

2 

2 

Gontinuous 

Yes 

10,272 

5. 1.1. 3 

Tamaki 

3 

3 

Gontinuous 

Yes 

145 

5. 1.1. 4 

Dias  PI 

30 

2 

Gontinuous 

Yes 

697 

5. 1.2.1 

Dias  r2 

30 

2 

Gontinuous 

No 

625 

5. 1.2.2 

Fonseca  FI 

2 

2 

Gontinuous 

Yes 

10,036 

5. 1.2.3 

Schaffer  F3 

1 

2 

Gontinuous 

Yes 

11,250 

5. 1.2.4 

Srinivas 

2 

2 

Gontinuous 

Yes 

697 

5. 1.2.5 

DTLZ7 

2 

2 

Gontinuous 

No 

36 

5. 1.2.6 

Disk  Brake 

4 

2 

Mixed 

Yes 

108 

5. 1.3.1 

The  algorithm  was  then  applied  to  an  airfoil  design  optimization  problem.  This 
problem  is  representative  of  real-world  multi-objective  aircraft  design  problems.  (Also 
see  [35],  [77],  [90],  and  [82].)  Additionally,  objectives  of  engineering  design  optimiza¬ 
tion  problems  are  often  either  subject  to  measurement  error  or  must  be  estimated 
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with  simulations.  Examples  of  simulation  used  in  aircraft  design  can  be  found 
in  [91],  [81],  [36],  and  [35].  Results  of  the  application  are  shown  in  Section  5.2. 


5. 1  Test  Results 

5.1.1  Continuous  Multi- Objective  Test  Problems. 


5. 1.1.1  Viennet  j.  As  a  proof-of- concept.  Implementation  1  of  the 
algorithm  (see  Section  4.5)  was  built  and  tested  on  a  problem  given  by  Viennet  and 
Marc  [149]: 


min  Fi{Xi,X2) 
-f2(Vi,  X2) 
A3(Xi,X2) 

subject  to 


(^1-2)^  ,  (X2  +  l)^  , 

(Xi  +  X2-3)^  {2X2 -X^f  _ 

1 75  1 7 

(3Xi-2X2  +  4)2  (Xi-X2  +  1)2 
- - h  15 


AXi  -j-  X2  —  4^0 


-ATi  -  1  <  0 


Xi-  X2 -  2  <0 


Xi,X2e[-4,+4]2. 


This  test  problem  was  modihed  by  adding  uniformly  distributed  random  noise 
between  0  and  1  to  each  objective  function  evaluation  to  simulate  the  algorithm’s 
use  on  a  stochastic  systemh  The  algorithm  was  tested  on  the  initial  problem  over  a 
range  of  aspiration  and  reservation  levels  using  three  different  experimental  designs: 
CCD  with  59  design  points,  BED  with  54  design  points,  and  FED  with  3  objectives 
set  at  4  levels  within  the  region  of  interest,  resulting  in  4,096  design  points.  Five 
replications  were  used  at  each  design  point.  Each  run  took  less  than  a  minute  (with 
500  function  evaluations). 

^This  is  true  of  all  further  test  problems  as  well  and  thus  will  not  be  mentioned  in  each  further 
test  problem. 
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(b)  Initial  Test  Results  for  Test  Problem  with  Added  Noise  using  Full  Factorial  Design 


Figure  5.1: 


Comparison  of  Initial  Test  Results  to  Published  Solution 
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Initial  results  are  shown  in  Figure  5.1.  The  initial  runs  fall  inside  the  pub¬ 
lished  Pareto  set  (see  Figure  5.1(a)),  implying  that  Stage  1  of  the  algorithm  is  indeed 
converging  to  Pareto  solutions.  Because  this  problem  contains  3  objective  functions, 
pairwise  examinations  of  the  experimental  Pareto  fronts  are  used.  These  appear  to 
approach  a  reasonable  approximation  to  the  actual  pairwise  Pareto  fronts.  However, 
the  effect  of  random  noise  appears  to  be  more  prominent  in  objective  2.  This  is  ex¬ 
pected  because  the  random  noise  was  not  scaled  and  is  much  larger  in  comparison  to 
values  of  objective  2  than  compared  to  those  of  objectives  1  and  3.  In  the  remaining 
test  problems,  noise  for  each  objective  was  scaled  to  approximately  one  percent  of 
its  maximum  objective  function  value  (he.,  value  at  the  nadir  point)  to  account  for 
differences  in  the  relative  sizes  of  the  objective  function  values,  so  that  each  is  affected 
similarly  by  the  noise. 

In  general,  as  the  magnitude  of  the  noise  grows,  the  number  of  iterations  re¬ 
quired  before  the  algorithm  terminates  may  become  prohibitively  large.  As  discussed 
in  Section  4.4.1,  one  of  the  stopping  criteria  in  the  subproblems  is  the  signal-to-noise 
ratio  squared  (AT  in  Equation  (4.1)).  Assuming  a  computational  budget  as  a  func¬ 
tion  of  this  value  AT,  the  sample  variance  should  be  less  than  ATtj^,  where  SI  is  the 
desired  indifference  zone  parameter.  Additionally,  because  the  number  of  iterations 
is  necessarily  finite,  the  filter  as  described  in  Section  3.2. 1.5  was  added  to  all  further 
implementations  (Implementations  2,  4,  6,  and  7  in  Table  4.1)  to  prevent  potentially 
dominated  points  from  entering  the  efficient  set. 
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5. 1.1. 2  Viennet  3  Test  Problem.  The  Viennet  3  test  problem  [149]  is 
linearly  constrained  and  mnlti-objective  and  is  given  by 

min  Fi{x,  y) 

F2{x,y) 

Fz{x,y) 

snbject  to 

—3  <  x,y  <3. 

Graphical  depictions  of  the  objective  fnnctions  can  be  seen  in  Fignre  5.2.  Using  im¬ 
plementation  2  (see  Section  4.5),  the  SMOMADS  algorithm  fonnd  an  approximation 
to  the  Pareto  optimal  set  and  front,  shown  in  Fignres  5.3(c)  and  5.3(d)  respectively, 
by  varying  the  aspiration  and  reservation  levels  nsing  a  fnll  factorial  design  with  3 
objectives  set  at  4  levels,  resnlting  in  4,096  design  points.  Becanse  solntions  retnrned 
by  the  snbproblem  may  be  hltered  and  not  enter  the  Pareto  set,  hve  replications  were 
nsed  at  each  design  point. 

When  compared  to  the  pnblished  resnlts  from  [148,149]  shown  in  Fignres  5.3(a) 
and  5.3(b),  it  is  clear  that  the  SMOMADS  algorithm  was  able  to  prodnce  a  good 
approximation  of  the  pnblished  solntion.  Althongh  each  snbproblem  replication  was 
relatively  fast,  the  large  nnmber  of  design  points  associated  with  the  fnll  factorial 
design  led  to  several  honrs  of  compntation.  This  time  conld  be  rednced  by  nsing  a 
central  composite  design  to  examine  the  entire  space  and  then  focns  in  areas  where 
gaps  may  exist  in  the  approximation  to  the  Pareto  front.  Additionally,  test  metrics, 
dehned  in  Section  3.2. 1.6,  were  nsed  to  determine  the  qnality  of  the  Pareto  set  after 
the  rnn  was  complete.  These  valnes  are  shown  in  Table  5.3.  However,  becanse  these 
metrics  are  generally  nsed  to  compare  Pareto  sets,  it  wonld  be  nsefnl  in  determining 
if  the  algorithm  shonld  terminate  before  the  entire  experimental  design  had  been 
execnted  if  a  desired  qnality  in  the  Pareto  set  approximation  or  threshold  valne  had 


=  0.5  (x^  -f-  y^)  +  sin  (x^  -f  y"^) 

_  {3x-2y  +  Af  ^  {x-y  +  lf  ^ 
— - - 1  - r 


{x^  +  y‘^  +  1) 
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Viennet  3  —  Objective  1  with  Projections 


Viennet  3  —  Objective  2  with  Projections 


(a)  Objective  Function  1  with  Projections 
onto  the  axes 


VariablB  2  ^  -4  Variable  1 

(b)  Objective  Function  2  with  Projections 
onto  the  axes 


(c)  Objective  Function  3  with  Projections  (d)  All  Objective  Functions 

onto  the  axes 


Figure  5.2:  Graphical  Depiction  of  the  Viennet  3  Test  Problem 
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(a)  Pareto  Optimal  Set  for  Deterministic 
Test  Problem,  Figure  6  in  [149] 


(b)  Pareto  Front  for  Deterministic  Test  Prob¬ 
lem,  Figure  C.36  in  [148] 


Viennet  3  Test  Problem  —  Experimental  Pareto  Front 


(c)  Experimental  Pareto  Optimal  Set 


(d)  Experimental  Pareto  Front 


Figure  5.3:  Comparison  of  Experimental  Results  to  Published  Solution  for  Viennet 
3  Test  Problem 
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been  obtained  (see  Section  4.4.2).  Fntnre  research  in  this  area  is  recommended  in 
Section  6. 2. 3. 2. 


Poloni  -  Objective  1  with  Projections 


Poioni  —  Objective  2  with  Projections 


(a)  Objective  Function  1  with  Projections 
onto  the  axes 


(b)  Objective  Function  2  with  Projections 
onto  the  axes 


Poloni  Test  Problem  —  Both  Objectives 


CM 

■a 

c 

cb 


S' 

O 


Variable  1 


(c)  Both  Objective  Functions 


Fignre  5.4:  Graphical  Depiction  of  the  Poloni  Test  Problem 


5. 1.1.3  Poloni 
Carlo  Poloni  in  1995  [117], 


Test  Problem.  The  Poloni  test  problem,  introdnced  by 
contains  two  nonlinear  objectives  and  simple  bonnds  on 
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the  variables.  It  is  given  by 


min  F{x,y)  =  {fi{x,y),  f2{x,y)) 
snbject  to 

— TT  <  x,y  <  n, 

where 

I/)  =  1  +  (^1  ~  BiY  +  {A2  —  B2Y 
f2ix,y)  =  (a;  +  3)2  +  (y  +  if 

Ai  =  0.5  sin(l)  —  2  cos(l)  +  sin(2)  —  1.5  cos(2) 

A2  =  1.5  sin(l)  —  cos(l)  +  2  sin(2)  —  0.5  cos(2) 

Bi  =  0.5  sin(a;)  —  2  cos(a;)  +  sin(|/)  —  1.5  cos(y) 

B2  =  1.5  sin(a;)  —  cos(a;)  +  2  sm{y)  —  0.5  cos{y), 

and  graphical  depictions  of  the  objective  fnnctions  are  shown  in  Fignre  5.4^.  Im¬ 
plementation  2  (see  Section  4.1)  of  the  SMOMADS  algorithm  was  nsed  to  hnd  the 
approximations  of  the  Pareto  optimal  set  and  Pareto  front  shown  in  Fignres  5.5(c) 
and  5.5(d),  respectively.  As  with  the  Viennet  3  test  problem,  the  objective  space 
was  divided  and  a  fnll  factorial  design  nsed.  In  this  case,  however.  Implementation  2 
tended  not  to  hnd  the  Pareto  points  that  heavily  favored  one  objective  or  the  other. 
That  is,  solntions  fonnd  by  the  algorithm  tended  to  he  in  the  “compromise  area”  of 
the  cnrve  shown  in  Fignre  5.5(d)  (area  near  the  middle  of  the  Pareto  front  where 
individnal  objectives  are  neither  extremely  high  nor  extremely  low).  After  adjnst- 
ing  the  range  of  the  aspiration  and  reservation  levels,  additional  rnns  resnlted  in  the 
generation  of  new  points  on  the  lower  right  side  of  the  cnrve. 

The  solntion  set  was  then  compared  to  the  one  pnblished  in  the  doctoral  dis¬ 
sertation  of  Van  Veldhnizen  [148].  Van  Veldhnizen’s  approximations  of  the  Pareto 

^Note  that  the  objective  functions  used  here  as  a  minimization  problem  are  the  negatives  of  the 
of  those  used  in  the  original  maximization  problem.  To  compare  with  published  results,  the  results 
in  Figures  5.5(c)  and  5.5(d)  are  converted  to  those  of  the  original  maximization  problem. 
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optimal  set  and  Pareto  front  for  the  deterministic  problem  (i.e.,  withont  noise  added), 
which  was  generated  by  an  evolntionary  algorithm,  are  shown  in  Fignres  5.5(a)  and 
5.5(b),  respectively.  As  can  be  clearly  seen  from  Fignre  5.5(b),  this  Pareto  set  ap¬ 
proximation  contains  dominated  points.  In  contrast,  the  Pareto  set  approximation 
fonnd  via  SMOMADS  contains  few  (if  any)  dominated  points  as  shown  in  Fignre 
5.5(d).  Thongh  the  Van  Veldhnizen  set  contains  more  points  in  the  approximation 
and  covers  more  of  the  objective  space,  it  shonld  be  noted  that  the  SMOMADS  al¬ 
gorithm  was  not  rnn  exhanstively.  Additional  aspiration-reservation  level  pairs  can 
be  rnn  in  order  to  determine  additional  portions  of  the  Pareto  frontier.  The  limiting 
factor  is  the  nnmber  of  design  points  in  the  master  problem  which  can  be  evalnated. 
Thns,  as  previonsly  mentioned,  efficient  methods  for  searching  the  objective  space  can 
be  of  paramonnt  importance. 

5. 1.1. 4  Tamaki  Test  Problem.  The  Tamaki  test  problem  [148]  contains 
three  linear  objectives,  one  nonlinear  constraint,  and  nonnegativity  constraints  on  the 
variables.  Its  formnlation  is  given  by 

max  F{x,  y,  z)  =  y,  z),  f2{x,  y,  z),  ^x,  y,  z)) 

snbject  to 

x"^  +  y"^  +  z"^  <  1 
x,y,z>0 

where 

fi{x,y,z)  =  X 
f2{x,y,z)  =  y 
f3{x,y,z)  =  z. 

As  shown  in  Fignre  5.6(a),  its  feasible  region  is  the  nonnegative  portion  of  a  sphere. 
Becanse  the  objective  is  to  maximize  each  variable,  the  trne  Pareto  front  is  the  snrface 
of  the  sphere  within  the  feasible  region  (see  Fignre  5.7(a)). 
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PotanI  Pa«toOpdmi(  Sdulons 
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POIOM  PmMO  FrdM 


(a)  Pareto  Optimal  Set  for  Deterministic  (b)  Pareto  Front  for  Deterministic  Test  Prob- 

Test  Problem,  Figure  C.17  in  [148]  lem,  Figure  C.18  in  [148] 


Poloni  Test  Problem  —  Experimental  Pareto  Optimal  Set 


Poloini  Test  Problem  —  Experimental  Pareto  Front 


-70  -60  -50  -40  -30  -20  -10 

Objecitve  1 


(c)  Experimental  Approximation  of  the 
Pareto  Optimal  Set,  Determined  Via  Full 
Factorial  Design  of  the  Region  of  Interest 


(d)  Experimental  Approximation  of  the 
Pareto  Front,  Determined  Via  Full  Factorial 
Design  of  the  Region  of  Interest  Showing  Few 
Dominated  Points 


Figure  5.5:  Comparison  of  Experimental  Results  to  Published  Solution  for  Poloni 
Test  Problem 
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Implementations  6  and  7  of  the  SMOMADS  algorithm  were  used  to  hnd  an 
approximation  of  the  Pareto  optimal  set/front  shown  in  Figures  5.7(c)  and  5.7(d). 
As  with  previous  test  problems,  the  algorithm  tended  toward  “compromise  ”  regions 
of  the  objective  space.  However,  this  is  due  to  the  choice  of  experimental  design  in 
the  master  problem,  which  can  be  changed  to  use  a  different  implementation  of  the 
algorithm,  and  is  thus,  not  considered  an  algorithmic  shortcoming. 

5.1.2  Continuous  Bi- Objective  Test  Problems. 

5. 1.2.1  Dias  FI  Test  Problem.  The  Dias  FI  test  problem  is  a  bound 
constrained  bi-objective  (one  linear  function,  one  nonlinear  function)  problem  with 
30  decision  variables.  Its  formulation  is  given  by 

mmF{X)  =  {f,(X),MX)) 
subject  to 

0<  ATi  <  1,  i  =  1,2,...,30 

where 


fi(X)  =  X, 


X  =  [X„...,X3o]. 


Implementation  2  of  the  SMOMADS  algorithm  was  used  to  determine  the  approxima¬ 
tion  of  the  Pareto  front  (shown  in  Figure  5.8(b)),  as  follows.  First,  implementation  2 
was  run  in  the  range  between  the  utopia  and  nadir  points  in  the  objective  space.  The 
graph  of  the  Pareto  front  was  observed  and  areas  with  gaps  and  under-representation 
were  determined.  Further  runs  (using  Implementation  4)  were  then  focused  on  those 
areas  by  setting  the  aspiration  and  reservation  range  to  a  small  area  around  the  miss¬ 
ing  regions.  When  compared  to  the  published  solution  to  the  deterministic  version 
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Tamato  ~  Feasible  Region 


Tamaki  —  Paraeto  Optimal  Set/Pareto  Front 


Variable/Objective  2 


Variable/Objective  1 


(a)  Feasible  Region 


(b)  Theoretical  Pareto  Front 


Figure  5.6:  Theoretical  Solutions  for  the  Tamaki  Test  Problem 


(a)  Pareto  Optimal  Set  for  Deterministic  (b)  Pareto  Front  for  Deterministic  Test  Prob- 

Test  Problem,  Figure  D.19  in  [148]  lem,  Figure  D.20  in  [148] 


Tamaki  Test  Problem  —  Pareto  Optimal  SehFront 


(c)  Experimental  Approximation  of  the 
Pareto  Optimal  Set  and  Front 


Tamaki  Test  Problem  —  Pareto  Optimal  Set/Front 


0.8- 


■|  0.6- 
O 


0.2- 


VariabiG/Objective  1  Variabla/ObJeMive  2 

(d)  Experimental  Approximation  of  the 
Pareto  Set  and  Eront 


Figure  5.7:  Comparison  of  Experimental  Results  to  Published  Solution  for  Tamaki 
Test  Problem 
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of  the  problem  found  in  [42],  the  approximation  determined  by  SMOMADS  is  quite 
good,  having  nearly  the  same  spread  and  fewer  points  that  are  actually  dominated. 


Dias  Gamma  1  Test  Problem  -  Experimental  Pareto  Front 


(a)  Pareto  Front  for  Deterministic  Test  Prob-  (b)  Experimental  Pareto  Front 

lem,  Figure  2  in  [42] 


Figure  5.8:  Comparison  of  Experimental  Results  to  Published  Solution  for  Dias’ 

ri  Test  Problem 


5. 1.2. 2  Dias  V2  Test  Problem.  The  Diaz  r2  problem  is  identical  to 
the  Diaz  PI  problem,  except  for  a  change  in  the  second  objective.  It  is  given  by 


minF(X)  =  (/i(X),/2(X)) 

subject  to 

0  <  ATi  <  1,  i  =  1,2,. . .  ,30 

where 

h(X)  =  X, 


f2{X)  = 


^  A  a:,- 

1  +  9E 


^^2  V(^-i) 


X  =  [X,,...,X3o]. 


/ 


MX) 


\ 


M 

1  +  9E 

i=2 


X,  \ 
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Gamma  2  Test  Problem  —  Experimental  Pareto  Front 


(a)  Pareto  Front  for  Deterministic  Test  Prob-  (b)  Experimental  Pareto  Front 

lem,  Figure  3  in  [42] 

Figure  5.9:  Comparison  of  Experimental  Results  to  Published  Solution  for  Dias’ 

T2  Test  Problem 

5. 1.2. 3  Fonseca  FI  Test  Problem.  The  Fonseca  FI  test  problem  con¬ 
tains  two  exponential  objectives  and  two  variables  with  upper  and  lower  bounds. 
Graphs  of  both  objective  functions  are  shown  in  Figures  5.10(a)  and  5.10(b),  and  the 
problem  formulation  is  given  by 

min  fi{Xi,  X2)  =  1-  exp  (-(Xi  -  1)^  -  (X2  +  1)^) 

MX,,  X2)  =  1  -  exp  (-(Xi  +  1)2  -  (X2  -  1)2) 

subject  to 

-2  <  Xi  <  2 
-2  <  X2  <  2. 


The  theoretical  Pareto  optimal  set  (as  shown  in  Figure  5.11(a))  can  be  easily  de¬ 
termined  to  be  the  line  between  points  (0,1)  and  (1,0)  [79].  The  corresponding 
theoretical  Pareto  front  is  shown  in  Figure  5.11(b). 

Implementations  2  and  4  of  the  SMOMADS  algorithm  were  used  to  determine 
the  experimental  Pareto  optimal  set  and  front  as  shown  in  Figures  5.11(c)  and  5.11(d), 
respectively.  As  shown  in  the  graphs,  the  algorithm  found  a  very  good  approximation 
in  the  “compromise  region”  as  well  as  the  extreme  solutions.  Gaps  do  exist  in  the 


approximation  that  could  be  examined  further  via  the  master  problem  termination 
algorithm  outlined  in  Section  4.4.2. 

5. 1.2. 4  Schaffer  F3  Test  Problem.  The  Schaffer  F3  test  problem  con¬ 
tains  one  quadratic  objective  function,  one  nonsmooth  piecewise  linear  objective  func¬ 
tion,  and  one  decision  variable  with  simple  bounds.  The  problem  formulation  is  given 
by 

min  F{X)  =  (fi(X),MX)) 
subject  to 

-8  <  X  <  8 

where 

-X,  X  <  1 

-2  +  X,  1  <  X  <  3 

/i(^)  =  <^ 

4-X,  3<X<4 

-4  -f  X,  4  <  X 

\ 

/2(X)  =  (X-5)2. 

Graphical  depictions  of  the  objectives  are  shown  in  Figures  5.12(a)  and  5.12(c).  The 
published  Pareto  optimal  set  is  also  shown  in  Figure  5.12(a).  The  true  Pareto  front 
(derived  from  the  published  true  Pareto  set)  is  shown  in  Figure  5.12(b). 

Implementations  2  and  4  were  used  to  determine  the  experimental  Pareto  op¬ 
timal  set  and  front  shown  in  Figures  5.12(c)  and  5.12(d),  respectively.  As  depicted 
in  the  graphs,  the  experimental  solution  nearly  perfectly  matches  the  published  the¬ 
oretical  solution. 

5. 1.2. 5  Srinivas  Test  Problem.  The  Srinivas  test  problem  contains  two 
nonlinear  objective  functions,  two  nonlinear  constraints,  and  two  decision  variables 
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Fonsecas  Test  Problem  -  Objective  Functions 


Fonsecas  Test  Problem  -  Objective  Functions 


Variable  t 


(a)  Objective  Functions 


(b)  Objective  Functions  -  Alternate  View 


Figure  5.10:  Graphical  Depiction  of  the  Fonseca  FI  Test  Problem 


Fonseca’s  FI  Test  Problem  —  Pareto  Optimal  Set 


Fonseca's  FI  Test  Problem  —  Pareto  Front 


(a)  Pareto  Optimal  Set  for  Deterministic  (b)  Pareto  Front  for  Deterministic  Test  Prob- 
Test  Problem  [79]  lem 


Fonsecas  FI  Test  Problem  -  Experimental  Pareto  Optimal  Set 


^  i  k  :  ;• 


■‘M: 


-1.5  -1  -0.5  0  0.5  1  1.5 

Variable  1 


(c)  Experimental  Pareto  Optimal  Set 


Fonsecas  FI  Test  Problem  -  Experimental  Pareto  Front 


0.2  0.4  0.6  0.8  1 

Objective  1 

(d)  Experimental  Pareto  Front 


Figure  5.11:  Comparison  of  Experimental  Results  to  Published  Solution  for  Fonseca 
FI  Test  Problem 
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Schaffer's  F3  Test  Problem — Pareto  Front 


(a)  Pareto  Optimal  Set  for  Deterministic 
Test  Problem,  Figure  4  in  [79] 


(b)  Pareto  Front  for  Deterministic  Test  Prob¬ 
lem 


(c)  Experimental  Pareto  Optimal  Set 


(d)  Experimental  Pareto  Front 


Figure  5.12:  Comparison  of  Experimental  Results  to  Published  Solution  for  Schaffer 
F3  Test  Problem 


with  upper  and  lower  bounds.  Its  formulation  is  given  by 


min  fi{Xi,X2)  =  1  -exp(-(Xi  -  1)^  -  {X2  +  1)^) 
f2{Xi,X2)  =  1  -exp(-(Xi  +  1)2  -  (X2  -  1)2) 


subject  to 


Xf  +  X2  -  225  <  0 
Xi  -  3X2  +  10  <  0 


-20  <  Xi  <  20 


-20  <  X2  <  20, 
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and  graphical  depictions  of  the  feasible  region  and  objective  fnnctions  are  shown  in 
Fignre  5.13.  The  Pareto  optimal  set  and  front,  as  pnblished  by  Van  Veldhnizen,  are 
shown  in  Fignres  5.14(a)  and  5.14(b),  respectively  [148]. 

Implementations  6  and  7  were  nsed  to  determine  the  experimental  Pareto  op¬ 
timal  set  and  front  shown  in  Fignres  5.14(c)  and  5.14(d),  respectively.  As  depicted 
in  the  graphs,  the  Pareto  optimal  set  and  front  fonnd  via  the  SMOMADS  algorithm 
match  qnite  well  to  the  pnblished  solntion.  In  fact,  the  pnblished  solntion,  fonnd  via 
an  evolntionary  algorithm,  clearly  contains  many  dominated  solntions,  whereas  the 
experimental  solntion  contains  very  few. 


Srinivas  —  Feasible  Region  Srinivas  —  Objective  1  with  Projections 


(a)  Feasible  Region  (b)  Objective  Function  1  with  Projections 

onto  the  axes 


300 
250 
200 
I  150 
°  100 
50 
0 


(c)  Objective  Function  2  with  Projections  (d)  Both  Objective  Functions 

onto  the  axes 


Fignre  5.13:  Graphical  Depictions  of  the  Srinivas  Test  Problem 
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(a)  Pareto  Optimal  Set  for  Deterministic 
Test  Problem,  Figure  D.17  in  [148] 


Srlnlvas  Pareto  Front 


(b)  Pareto  Front  for  Deterministic  Test  Prob¬ 
lem,  Figure  D.18  in  [148] 
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Srinivas  Test  Problem  -  Experimental  Pareto  Optimal  Set 
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(c)  Experimental  Pareto  Optimal  Set 
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Srinivas  Test  Problem  -  Experimental  Pareto  Front 
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Objective  1 

(d)  Experimental  Pareto  Front 


Figure  5.14:  Comparison  of  Experimental  Results  to  Published  Solution  for  Srinivas 
Test  Problem 


5. 1.2. 6  DTLZ7  (2  Objective)  Test  Problem.  The  DTLZ7  test  problem 
can  be  found  in  [119]  and  [41],  It  is  bi-objective  with  one  linear  and  one  nonlinear 
objective  function  and  two  decision  variables  with  lower  and  upper  bounds.  The 
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problem  formulation  is 


min  f,{X^,X2)  =  X, 


/2(Xi,X2)  =  (1  +  10X2) 


Xi 


1  +  lOX 


Xi  sin(87rXi) 

1  +  10X2 


subject  to 


0  <  Xi,  X2  <  1. 


The  pubished  Pareto  optimal  front,  as  determined  via  a  particle  swarm  algorithm 
[119],  is  shown  in  Figure  5.15(a).  Implementations  2  and  4  were  used  to  determine 
the  experimental  solution  shown  in  Figure  5.15(b).  When  compared  to  the  published 
solution,  the  experimental  solution  matches  quite  well. 


Swarm  Algorithm 


2  Objective  DTL7  Problem  —  Experimental  Pareto  Front 


0.8 
I  0.6 
I  0.4 
’  0.2 
0 

-0.2 

-0.4 

-0.6 


0.4  0.5 

Objeotive  1 


(a)  Pareto  Front  for  Deterministic  Test  Prob-  (b)  Experimental  Approximation  of  the 
lem,  Figure  1  in  [119]  Pareto  Front 


Figure  5.15:  Comparison  of  Experimental  Results  to  Published  Solution  for  2  Ob¬ 
jective  DTLZ7  Test  Problem 


5.1.3  Mixed  Variable  Bi-Objective  Test  Problems. 

5. 1.3.1  Disk  Brake  Design  Test  Problem.  This  mixed  variable  test 
problem  is  found  in  [119]  and  is  an  engineering  design  problem  for  a  disk  brake.  It 
contains  two  nonlinear  objective  functions,  two  linear  constraints,  three  continuous 
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variables  and  one  categorical  variable,  and  is  given  by 


max 


h{X,,X,,X3,X,)  =  (4.9  X  10-5)(X| 

19.82  X 

f2  Xi,X2,X3,X4 


..  2-x!){x, 

(9.82  X  10®)(X|  -X2) 


1) 


snbject  to 


0.4 


{X2  -  Xi)  -  20  >  0 

30  -  2.5(X4  +  1)  >  0 
X3 


1  - 


3.14(X|-X2) 
(2.22  X  10~3)X3(X|  -Xf) 


(X|  -  Xf) 
(2.66  X  10~^)X3X4{Xl  -  Xf) 
(^1  - 


>  0 


>  0 


-  900  >  0 


55  <  Xi  <  80 
75  <  X2  <  110 
1000  <  X3  <  3000 
2  <  X4  <  20 
X4  e  N. 


The  lone  categorical  variable  X4  represents  the  nnmber  of  disks  in  the  brake.  Thongh 
it  can  take  on  only  19  distinct  valnes,  meaning  that  exhanstive  ennmeration  is  not 
prohibitively  expensive,  the  point  of  this  test  problem  was  to  exercise  the  mixed 
variable  logic  in  the  solver.  Implementation  7  of  the  SMOMADS  algorithm  was 
nsed  to  determine  the  Pareto  front  shown  in  Fignre  5.16(b).  When  compared  to 
the  pnblished  solntion  (generated  via  a  particle  swarm  algorithm)  [119],  the  solntion 
generated  by  SMOMADS  is  nearly  identical. 


5.2  Engineering  Design  Optimization  Applieation 

After  the  algorithm  snccessfnlly  solved  the  pnblished  test  problems,  it  was  ap¬ 
plied  to  a  small  airfoil  engineering  design  optimization  problem.  Mathematically,  the 
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35 


Swarm  Algorithm 


‘^tbo 


CQjc 


1.5 

Function  1 


(a)  Pareto  Front  for  Deterministic  Test  Prob¬ 
lem,  Figure  6  in  [119] 


Disk  Brake  Test  Problem  —  Experimental  Pareto  Front 


0  0.5  1  1.5  2  2.5  3 

Objective  1 

(b)  Experimental  Approximation  of  the 
Pareto  Front 


Figure  5.16:  Comparison  of  Experimental  Resnlts  to  Pnblished  Solntion  for  Disk 
Brake  Test  Problem 


problem  is: 


^min^  [E  [q]  ,  (q)  ,  (q)  ,  -E  [q])  (5.12a) 

where 

Q  =  /i(Xi,X2,X3) 

C,  =  /2(Xi,X2,X3) 

snbject  to 

E  [ci]  >  0.612  (q)  <  0.0546 

E  [cd]  <  0.0172  (cd)  <  0.0061 

1^  [cm]  I  <  0.0888  (c^)  <  0.00827 

0  <  ATi  <  10  0  <  X2  <  10  12  <  273  <  100 

This  problem  was  adapted  from  one  originally  solved  by  Ciprian  et  al.  [29]  nsing  mnlti- 
criteria  decision  making  methods  and  game  theory.  The  decision  variables  Xi  and  X2 
represent  the  hrst  and  second  digits  in  the  NACA^  airfoil  nnmber  and  X3  represents 

^NACA  is  the  National  Advisory  Committee  for  Aeronautics,  established  by  an  act  of  the  United 
States  Congress  on  3  March  1915  [7].  NACA  ceased  to  exist  on  1  October  1958  and  was  replaced  by 
the  National  Aeronautics  and  Space  Administration  (NASA)  [108];  however,  the  NACA  numbering 


(5.12b) 


(5.12c) 
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the  third  and  fourth  digits.  Following  the  standard  NACA  airfoil  numbering  system, 
this  implies  that  Xi,  X2,  and  X3  represent,  respectively,  the  maximum  airfoil  camber 
(in  percentage  of  its  chord,  the  distance  between  the  leading  and  trailing  edges  of  an 
airfoil  measured  in  the  direction  of  normal  airflow  [156]),  distance  of  maximum  camber 
from  the  airfoil  leading  edge  (in  tens  of  percents  of  its  chord),  and  the  maximum 
thickness  of  the  airfoil  (in  percent  of  its  chord)  [157].  As  shown  in  Equation  (5.12a), 
the  objectives  are  to  minimize  the  expected  value  and  variance  of  the  drag  coefficient 
and  the  variance  of  the  lift  coefficient  and  to  maximize^  the  expected  value  of  the  lift 
coefficient  and  the  lift  and  drag  coefficients  are  assumed  to  be  black  box  functions  of 
the  decision  variables.  In  the  original  problem,  the  coefficients  of  lift  and  drag  are 
determined  using  the  Navier-Stokes  version  of  the  flow  solver  MUFLO  and  AIRFOIL 
(using  different  shape  coordinates  as  decision  variables).  Neither  of  these  solvers 
were  available  open  source,  so  the  problem  was  modified  slightly  so  that  it  could  be 
used  with  the  PABLO  solver  (which  stands  for  Potential  Flow  around  Airfoils  with 
Boundary  Layer  Coupled  One-way) ,  a  pedagogical  low-speed  airfoil  analysis  program 
written  in  MATLAB®  [152].  The  changes  made  to  the  design  problem  are  as  follows: 


1.  The  original  problem  used  two  10-degree  Bezier  curves  to  represent  the  airfoil 
shape  and  used  their  coordinates  as  decision  variables.  The  modified  problem 
used  the  NACA  airfoil  numbering  system  to  represent  the  airfoil  shape  and  used 
the  first,  second,  and  third/fourth  digits  as  decision  variables. 

2.  The  flow  solver  in  the  original  problem  used  Mach  number  as  an  input.  The 
solver  in  the  modified  version  uses  Reynolds  number  as  an  input.  However, 
since  the  Reynolds  number  is  directly  proportional  to  velocity  by 


^  7IT7 
1? 


pVsL  _  VsL 
fX  V 


(5.13) 


system  continues  to  be  used. 

^The  expected  value  of  the  lift  coefficient  is  maximized  by  minimizing  its  negative. 
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where 


Vs  -  mean  fluid  velocity 
L  -  characteristic  length 
yU  -  (absolute)  dynamic  fluid  viscosity 
p  -  fluid  density 

V  -  kinematic  fluid  viscosity:  v  =  p/p 


and 


a 


V 


(5.14) 


where 

M  -  Mach  number 
V  -  velocity 

7  -  ratio  of  specific  heats  (for  air  7  =  1.4) 

R  -  specific  gas  constant 
T  -  temperature, 

and  for  a  fixed  temperature  and  pressure,  velocity  can  be  determined  from  Mach 
number  via  Equation  (5.14)  [7,158].  Unfortunately,  the  temperature,  pressure, 
and  characteristic  length  of  the  airfoil  are  not  given  in  [29];  thus,  standard 
temperature  and  pressure  at  sea  level  and  the  characteristic  length  of  ten  feet 
are  assumed. 

3.  The  original  problem  was  near  the  transonic  range,  the  range  of  speeds  just 
below  and  above  the  speed  of  sound,  approximately  Mach  0.8  -  1.2.  The  flow 
solver  in  the  modified  problem  is  subsonic,  and  the  speeds  in  original  problem 
(mach  0.73  ±0.05)  approach  the  limits  of  a  low-speed  solver,  so  the  speed  range 
in  the  modified  problem  was  reduced  to  mach  0.53  ±  0.05. 
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In  both  the  original  and  modihed  problems,  the  objectives  are  to  be  minimized 
over  a  range  of  Mach  nnmbers  (Reynolds  nnmbers)  and  angle  of  attack  (2  ±  0.5)  in 
order  to  hnd  a  robnst  design.  (However,  in  the  modihed  problem  nsing  the  low-speed 
solver,  the  objectives  do  not  vary  across  the  range  of  Mach  nnmbers,  so  the  variation 
is  dne  only  to  angle  of  attack.)  The  decision  space  is  constrained  snch  that  the  new 
airfoils  can  perform  no  worse  than  the  RAE2822  airfoiR  and  the  minimnm  airfoil 
thickness  is  12%  of  the  chord.  (The  performance  valnes  for  the  RAE2822  airfoil  were 
taken  from  [29]  and  are  given  in  Eqnation  (5.12c)).  This  problem  has  been  snfficiently 
modihed  snch  that  direct  comparison  to  the  pnblished  solntion  is  no  longer  germane 
and  is  not  shown. 

Implementation  7  of  the  SMOMADS  algorithm  was  nsed  to  determine  the  ex¬ 
perimental  Pareto  front,  shown  in  Fignre  5.17.  In  Fignres  5.17(a)  and  5.17(b),  the 
feasible  region  is  enclosed  (below  and  to  the  left)  by  the  solid  lines. 

Becanse  the  problem  had  to  be  modihed  to  rnn  inside  the  MATLAB®  environ¬ 
ment,  comparisons  to  the  pnblished  solntion  wonld  be  inappropriate.  However,  some 
observations  regarding  the  experimental  solntion  follow.  The  SMOMADS  algorithm 
is  able  to  hnd  many  non-dominated  solntions,  bnt  each  rnn  took  a  signihcant  amonnt 
of  time  (approximately  3  days  for  100  design  points).  This  nnderscores  the  necessity 
of  a  parsimonions  experimental  design  in  the  master  problem  and  highlights  the  need 
for  follow-on  research  in  this  area. 

As  with  the  test  problems,  the  SMOMADS  algorithm  tended  to  hnd  solntions  in 
the  “compromise”  area  of  the  Pareto  front.  This  particnlar  type  of  problem  highlights 
the  notion  that  snch  solntions  are  particnlar  nsefnl  (see  Section  4.4.2).  For  example, 
a  Pareto  optimal  solntion  corresponding  to  high  lift,  bnt  also  high  drag,  may  resnlt 
in  an  aircraft  that  nses  excessive  amonnts  of  fnel.  Similarly,  one  with  very  low  drag, 
bnt  little  lift,  conld  resnlt  in  an  aircraft  that  cannot  carry  mnch  weight.  Thns, 
the  solntions  fonnd  by  the  algorithm  wonld  be  those  more  likely  to  be  favored  by  a 

^RAE  stands  for  the  Royal  Aircraft  Establishment  and  is  the  British  version  of  the  United  States’ 
NACA  airfoil  system  [7]. 
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decision  maker. 


Airfoil  Design  Problem  —  Experimental  Pareto  Front  (Objectives  1  and  4) 


0.02  r 
0.018  - 


RAE2822 


0.016  - 
0.014  - 
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0.004  - 
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(a)  Experimental  Pareto  Front  for  Objectives  1  and  4 

Airfoil  Design  Problem  —  Experimental  Pareto  Front 
(Objectives  2  and  3) 


-0.5 


0.03  0.04  0.05  0.06 

Var(CI) 

(b)  Experimental  Pareto  Front  for  Objectives  2  and  3 

Figure  5.17:  Experimental  Results  of  the  Airfoil  Engineering  Design  Problem 
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5. 3  Quality  Metrics 


The  Pareto  set  quality  metrics  (Pareto  spread,  number  of  distinct  choices,  and 
cluster)  discussed  in  Section  3.2. 1.6  were  calculated  for  all  the  test  problems.  The 
parameters  of  the  test  problems  necessary  for  the  calculations  are  given  in  Table 
5.2  and  the  values  of  the  test  metrics  are  given  in  Table  5.3.  Since,  as  previously 
mentioned,  these  metrics  could  be  useful  in  determining  when  to  terminate  the  master 
problem,  the  experimental  Pareto  set  for  each  test  problem  was  reduced  by  90%  by 
randomly  selecting  10%  of  the  Pareto  optimal  points  and  then  reevaluating  the  metrics 
in  order  to  demonstrate  how  they  might  be  used  with  threshold  values.  The  values 
of  the  metrics  on  the  reduced  sets  are  shown  in  Table  5.4,  where  it  is  evident  that  the 
test  metrics  do  indeed  improve  as  the  number  of  points  in  the  Pareto  set  increases. 


Table  5.2:  Parameters  Required  for  the  Quality  Metrics 


Problem 

n 

/I 

/I 

/I 

/f 

/I 

fi 

/I 

np 

.1  X  np 

Viennet  3 

10 

18 

0.2 

1 

15 

-0.2 

3467 

346 

Viennet  4 

7.5 

-11 

26 

3.3 

13 

15 

58990 

5899 

Poloni 

30 

50 

0 

0 

35178 

3417 

Tamaki 

0 

0 

0 

-1 

-1 

-1 

404 

41 

Dias  PI 

1 

1 

0 

0 

772 

77 

Dias  r2 

1.1 

1.1 

0 

0 

330 

33 

Fonseca  FI 

1.01 

1.01 

0 

0 

1066 

107 

Schaffer  F3 

1 

16 

-1 

0 

1958 

196 

Srinivas 

250 

10 

0 

-250 

996 

97 

DTLZ7 

0.85 

1.4 

0 

-0.6 

377 

38 

Disk  Brake 

2.75 

33 

0 

0 

289 

29 

Airfoil 

0.0115 

0.004 

0.00025 

-0.612 

0 

0 

0 

-2.5 

231 

24 

Because  the  sets  were  reduced  in  a  random  fashion,  the  quality  of  the  sets  of  the 
different  test  problems  were  not  reduced  equally,  meaning  that  improvement  is  more 
significant  in  some  cases  than  in  others.  This  also  provides  some  insight  into  the 
utility  of  the  metrics.  If  the  quality  metrics  are  tracked  after  each  design  point,  they 
could,  in  fact,  be  used  as  stopping  criteria,  because  relatively  few  (quality)  points  can 
produce  good  values  in  the  test  metrics.  For  example,  in  the  Viennet  3  and  Viennet 
4  test  problems,  values  for  Pareto  spread  corresponding  to  the  the  reduced  Pareto 
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optimal  set  are  nearly  as  high  as  those  for  the  full  set.  If  this  had  been  the  hrst 
points  generated,  the  algorithm  could  have  been  terminated  after  only  10%  of  the 
experimental  design  had  been  evaluated.  Further  research  in  this  area  is  suggested 
in  Section  6. 2. 3. 2. 


Table  5.3:  Quality  Metrics  Applied  to  Test  Problems 


Problem 

OS{P) 

OSi{P) 

OS2{P) 

OSsiP) 

OS4{P) 

NDC^iP) 

CL^iP) 

Viennet  3 

0.4496 

0.9072 

0.7010 

0.7070 

96 

36.1 

Viennet  4 

0.5524 

0.9797 

0.5750 

0.9806 

173 

344.8 

Poloni 

0.2054 

0.3660 

0.5613 

53 

666.7 

Tamaki 

0.5304 

0.8287 

0.8535 

0.7499 

344 

1.2 

Dias  PI 

0.6987 

0.8013 

0.8719 

125 

6.2 

Dias  r2 

0.8619 

0.9182 

0.9388 

61 

5.4 

Fonseca  FI 

0.9942 

0.9975 

0.9968 

150 

7.1 

Schaffer  F3 

0.9933 

0.9936 

0.9997 

98 

20.0 

Srinivas 

0.7614 

0.8929 

0.8528 

155 

6.4 

DTLZ7 

0.8668 

0.9675 

0.8959 

76 

5.0 

Disk  Brake 

0.3830 

0.8294 

0.4618 

87 

3.3 

Airfoil 

0.1528 

0.4415 

0.4462 

0.8911 

0.8703 

43 

5.4 

Table  5.4: 

Problem 

Quality  Metrics  Applied  to  Reduced  Solution  Sets  of 
OS{P)  OSiiP)  OS2{P)  OSaiP)  OS^P)  NDC^{P) 

the  Test  Problems 
CL^iP) 

Viennet  3 

0.4296 

0.8977 

0.6790 

0.7048 

42 

8.7 

Viennet  4 

0.4245 

0.9303 

0.5496 

0.8302 

136 

43.3 

Poloni 

0.0164 

0.2925 

0.0559 

36 

98.0 

Tamaki 

0.0211 

0.3105 

0.3025 

0.2249 

40 

1.0 

Dias  FI 

0.6908 

0.7951 

0.8689 

54 

1.4 

Dias  r2 

0.8132 

0.8831 

0.9209 

19 

1.7 

Fonseca  FI 

0.9927 

0.9975 

0.9952 

57 

1.9 

Schaffer  F3 

0.9796 

0.9803 

0.9992 

69 

2.8 

Srinivas 

0.5817 

0.7668 

0.7586 

46 

2.1 

DTLZ7 

0.6557 

0.9396 

0.6978 

20 

1.9 

Disk  Brake 

0.3830 

0.8294 

0.4618 

26 

1.1 

Airfoil 

0.0175 

0.3671 

0.4462 

0.8897 

0.1203 

14 

1.7 

5.4  Efficiency  of  the  Algorithm 

The  computational  efficiency  of  SMOMADS  is  dependent  on  several  factors, 
which  are  determined  by  choices  made  during  implementation  of  the  algorithm.  Dis- 
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cussion  of  these  factors  and  observations  from  the  test  problems  follow. 

5.4- 1  Experimental  Design.  In  the  master  problem,  the  choice  of  exper¬ 
imental  design  directly  affects  the  compntation  time  reqnired  by  algorithm.  For 
example,  in  the  hrst  test  problem  (Viennet  4),  the  FFD  was  chosen,  resnlting  in  4,096 
individnal  snbproblems,  with  5  replications.  Compntational  time  for  each  replication 
took  considerably  less  than  a  minnte,  bnt  total  compntational  time  was  several  honrs. 
For  problems  with  expensive  or  lengthy  fnnction  evalnations,  snch  a  design  wonld  be 
prohibitively  expensive.  This  was  highlighted  by  snbseqnent  test  problems.  For 
example,  the  FFD  chosen  for  the  Poloni  test  problem  resnlted  in  10,000  individnal 
snbproblems  with  5  replications.  The  compntational  time  reqnired  was  over  72  honrs. 
If  a  CCD  had  been  chosen  instead,  the  nnmber  of  snbproblems  reqnired  wonld  have 
been  rednced  by  3  orders  of  magnitnde;  therefore,  this  design  was  nsed  in  later  test 
problems  in  a  seqnential  experimentation  scheme  (see  Section  4.1).  The  nnmber  of 
test  points  for  each  problem  and  whether  each  problem  was  solved  nsing  seqnential 
experimentation  are  shown  in  colnmns  6  and  5  of  Table  5.1,  respectively.  Every 
test  point  was  replicated  5  times,  except  for  the  hrst  CCD  (36  test  points)  of  the  disk 
brake  problem  that  was  replicated  15  times.  Farther  research  regarding  experimental 
designs  is  snggested  in  Section  6.2.3. 

5. 4-  2  Termination  Criteria.  The  overall  compntational  time  reqnired  by 
SMOMADS  is  heavily  dependent  on  the  compntational  efficiency  of  the  snbprob¬ 
lems®.  The  compntation  time  reqnired  by  the  snbproblems  depends  on  algorithmic 
parameters  of  the  two  individnal  solvers.  These  parameters  were  set  to  snggested 
defanlt  valnes  [3] ,  with  the  exception  of  the  maximnm  nnmber  of  iterations  allowed  in 
the  snbproblem,  which  was  set  to  500.  Limiting  the  nnmber  of  iterations  increased 
the  compntational  efficiency  by  stopping  snbproblems  that  might  be  converging  at 
a  slower  rate  and  thereby  snbstantially  increasing  the  overall  compntation  time.  In 

®The  number  of  objective  function  evaluations  required  by  the  solvers  is  not  captured  with  current 
code  and  could  not  be  reported.  Modifications  to  the  code  are  suggested  in  Section  6. 2.4.1 
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addition  to  directly  impacting  compntation  time,  limiting  iterations  also  indirectly 
impacts  compntational  efficiency  as  described  in  the  next  section. 

5.4.3  Solution  Quality  and  Filtered  Points.  In  the  limit,  the  snbproblems 
always  converge  to  Pareto  stationary  points;  however,  as  discnssed  in  Section  3. 2. 1.5, 
with  a  hnite  nnmber  of  iterations,  the  solver  may  terminate  at  a  dominated  point. 
Thns,  if  a  maximum  nnmber  of  iterations  is  stipnlated,  it  increases  the  likelihood 
of  prematnre  termination  and  directly  affects  solntion  qnality  in  the  snbproblems. 
Dominated  points  are  hltered  and  not  allowed  to  enter  the  efficient  set,  so  the  hnal 
solntion  qnality  is  not  affected;  however,  compntational  effort  is  wasted  for  each  point 
that  mnst  be  discarded.  Therefore,  the  maximnm  nnmber  of  iterations  mnst  be  set 
low  enongh  that  dnbions  snbproblems  are  abandoned,  yet  not  so  low  that  snbproblems 
are  sqnandered.  Unfortnnately,  this  can  be  problem-specihc.  For  example,  in  the 
Viennet  3  test  problem,  83%  of  points  were  hltered,  bnt  in  the  Poloni  test  problem, 
only  30%  of  solntions  did  not  enter  the  efficient  set. 

5. 5  Summary 

In  this  chapter,  the  accnracy  of  the  SMOMADS  algorithm  was  verihed  by  its 
nse  on  a  set  of  test  problems  with  known  solntions  and  an  airfoil  design  problem. 
The  algorithm  generally  prodnces  nnmerical  resnlts  qnite  similar  to  those  in  the  lit- 
eratnre.  In  the  design  problem,  it  was  able  to  generate  many  Pareto  optimal  points. 
Additionally,  factors  affecting  compntational  efficiency  were  discnssed  and  applicable 
examples  from  the  test  problems  were  illnstrated. 
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VI.  Conclusions  and  Recommendations 


A  research  approach  extending  stochastic  and  mnlti-objective  solntion  methods  to 
one  applicable  to  optimization  problems  having  both  characteristics  has  been  pre¬ 
sented.  The  new  method  fnrther  extends  the  generalized  pattern  search  methods  that 
previonsly  existed  for  single-objective  optimization  of  mixed- variable  systems  having 
stochastic  responses  to  the  mnlti-objective  case.  Snch  problems  are  typically  enconn- 
tered  when  one  desires  to  optimize  systems  with  mnltiple,  often  competing,  objectives 
that  do  not  have  a  closed  form  representation  and  mnst  be  estimated  via  simnlation. 
A  two-stage  method  was  presented  that  combines  generalized  pattern  search/rank¬ 
ing  and  selection  (GPS-R&S)  developed  for  single-objective  stochastic  problems  and 
Mixed  Variable  Mesh  Adaptive  Direct  Search  (MVMADS)  with  three  mnlti-objective 
methods:  interactive  techniqnes  for  the  specihcation  of  aspiration/reservation  levels, 
scalarization  fnnctions,  and  mnlti-objective  ranking  and  selection. 

Convergence  analysis  was  condncted  showing  that  in  each  snbproblem,  a  snbse- 
qnence  of  iterates  converges  with  probability  one  to  a  stationary  point  appropriately 
dehned  in  the  mixed-variable  domain,  and  that  if  that  solntion  is  also  optimal,  then 
it  is  Pareto  optimal.  Algorithmic  implementations  of  the  hrst  stage  have  been  de¬ 
veloped  and  compntational  experiments  condncted  on  standard  mnlti-objective  test 
problems  with  random  noise  added  to  the  objective  fnnction  evalnations  to  simnlate 
measnrement  error  or  the  nse  of  a  simnlation.  Finally,  the  algorithm  was  applied 
to  an  aeronantical  engineering  design  optimization  problem.  The  original  contribn- 
tions  of  this  dissertation  research  are  snmmarized  in  Section  6.1  and  fntnre  research 
directions  are  proposed  in  Section  6.2. 

6. 1  Contributions 

This  research  provided  the  following  original  contribntions  to  the  held  of  Oper¬ 
ations  Research. 
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6.1.1  General  Solution  Method  for  Multi- Objective  Optimization  of  Mixed  Vari¬ 
able  Systems  having  Stochastic  Responses.  The  primary  contribution  of  this  re¬ 
search  is  the  development  of  the  hrst  provably  convergent  algorithm  for  solving  multi¬ 
objective,  stochastic  optimization  problems  over  a  mixed-variable  domain.  Though 
other  methods  exist  for  problems  having  one  or  two  of  these  characteristics,  no  existing 
method  was  applicable  to  those  problems  exhibiting  all  three. 

6.1.2  Gonvergence  Analysis.  The  second  contribution  of  this  research  is 
a  rigorous  convergence  analysis  of  the  algorithm.  The  new  algorithm  was  shown 
to  converge  almost  surely  to  stationary  points  appropriately  dehned  in  the  mixed- 
variable  domain.  It  was  further  shown  that  if  such  points  are  optimal,  they  are 
also  Pareto  optimal.  In  proving  convergence,  it  was  necessary  to  hll  a  void  in  the 
existing  convergence  analysis  for  MADS,  i.e.,  to  show  that  a  mixed- variable  version 
of  MADS  converges  almost  surely  to  stationary  points  appropriately  dehned  in  the 
mixed  variable  case. 

6.1.3  Algorithmic  Implementation  Developed  and  Testing.  The  third  contri¬ 
bution  of  this  research  was  the  development  of  seven  implementations  of  the  algorithm 
that  were  tested  on  multi-objective  optimization  problems  from  the  literature.  These 
problems  had  random  noise  introduced  into  their  objective  function  evaluations  to 
simulate  random  responses.  Even  in  the  presence  of  noise,  the  algorithm  was  shown 
to  produce  results  comparable  to  the  published  solutions  that  were  generated  with 
deterministic  objective  function  evaluations. 

6.1.4  Application  Heuristic  Stopping  Griteria.  Heuristic  stopping  criteria 
were  developed  for  the  portion  of  the  algorithm  referred  to  as  the  master  problem. 
An  iterative  process  was  presented  (see  Section  4.4.2)  for  determining  when  an  ade¬ 
quate  approximation  to  the  Pareto  set  had  been  found.  Metrics  were  suggested  to 
substantiate  the  quality  of  the  set  and  for  termination. 
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6.2  Future  Research 


6.2.1  Convergence  Results  for  MVMADS-RS.  Convergence  resnlts  for  the 
stochastic,  mnlti-objective,  nonlinearly  constrained  case  have  not  been  rigoronsly 
proven  and  depend  on  Conjectnre  3.3.10.  A  formal  proof  of  this  conjectnre  is  re- 
qnired  for  completeness. 

6.2.2  Stage  2.  As  discnssed  in  Section  3.2.2,  some  efficient  points  may  not 
be  fonnd  via  interactive  specihcation  of  aspiration / reservation  levels  and  scalarization 
functions  of  the  type  discnssed  in  Section  3.1.4.  Thns,  a  second  (optional)  stage  was 
snggested  for  those  cases  in  which  missing  points  may  pose  a  particnlar  problem.  In 
this  stage,  a  discrete  mesh  (similar  to  that  nsed  for  MGPS-RS  and  MADS)  aronnd 
the  cnrrent  efficient  points  wonld  be  determined  and  then  a  mnlti-objective  ranking 
and  selection  algorithm  wonld  be  nsed  to  check  for  new  efficient  points  on  the  mesh. 
In  this  research,  the  algorithm  performed  so  well  on  the  test  problems  withont  the 
second  stage,  that  adding  it  was  not  necessary;  however,  fntnre  research  is  snggested 
for  those  cases  in  which  it  may  be  reqnired. 

6.2.3  Search  Methods  in  the  Objective  Space.  Thongh  the  proposed  method 
is  applicable  to  problems  with  any  nnmber  of  objectives,  in  practice,  the  dimension¬ 
ality  in  the  experimental  design  limits  searching  in  the  objective  space.  Thns,  an 
efficient  manner  of  searching  the  objective  space  is  needed,  and  the  following  possible 
research  areas  are  snggested. 

6.2.3. 1  Response  Surface  Methodology.  As  snggested  in  Section  4.1, 
the  problems  stndied  here  can  be  thonght  of  as  approximating  a  response  snrface 
(the  Pareto  front)  with  aspiration  and  reservation  levels  as  the  decision  variables 
nsing  experimental  design  methods.  In  this  research,  only  three  experimental  designs 
were  stndied:  fnll  factorial  design,  central  composite  design,  and  Box-Behnken  design. 
Other  methods  and  principles  from  experimental  design  may  be  as  nsefnl  in  efficiently 
approximating  the  Pareto  front.  Thns  a  more  in  depth  stndy  of  these  methods  is 
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suggested. 


6. 2. 3. 2  Quality  Metrics  for  the  Pareto  Set.  Currently,  the  experimen¬ 
tal  design  in  the  master  problem  of  the  solution  algorithm  is  run  exhaustively.  In 
some  cases  though,  an  adequate  Pareto  front  approximation  may  be  obtained  before 
the  entire  design  is  completed.  Thus,  if  the  quality  of  the  experimental  Pareto  set 
could  be  evaluated  at  various  points  and  the  search  terminated  at  a  predetermined 
threshold,  the  efficiency  of  the  search  method  could  be  improved. 

6. 2. 3. 3  Improved  Filter.  The  hlter  applied  in  this  research  used  a 
simple  comparison  of  objective  function  values  to  determine  if  a  point  was  dominated. 
Because  of  the  stochastic  nature  of  the  problem,  this  could  potentially  lead  to  some 
dominated  points  entering  the  set  or  non-dominated  points  being  unduly  excluded. 
The  experimental  results  in  Section  5.1  indicated  that  error  in  the  efficient  set  when 
compared  to  known  results  was  insignihcant;  however,  it  is  possible  that  in  other 
optimization  problems,  such  a  hlter  would  not  be  adequate.  One  suggestion  is  to 
add  a  multi-objective  ranking  and  selection  [85]  algorithm,  in  place  of  the  simple 
comparison  currently  in  the  hlter  to  determine  if  a  point  is  dominated  (see  Sections 
3.2.2. 1  and  3.1.5). 

Additionally,  the  hlter  only  determines  if  a  point  is  dominated  before  it  enters 
the  efficient  set.  Once  in  the  set,  it  is  never  again  considered.  It  is  possible  that  an 
erroneous  point  could  enter  the  efficient  set  early  in  the  algorithm’s  run.  Although  the 
experimental  results  in  this  study  did  not  indicate  signihcant  problems  with  solution 
accuracy,  it  is  recommended  that  the  hlter  also  determine  if  an  entering  efficient  point 
dominates  any  existing  efficient  solutions. 

6. 2. 3.4  Master  Problem  Termination  Criteria.  Though  heuristic  stop¬ 
ping  criteria  have  been  developed  here  (see  Section  4.4.2),  they  could  be  rehned  based 
on  response  surface  methodology  (Section  6.2.3)  and  quality  metrics  (Section  6. 2. 3. 2). 
With  an  astute  experimental  design,  a  more  thrifty  search  of  the  objective  space  could 
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be  conducted.  Furthermore,  with  insightful  quality  metrics  applied  after  each  design 
point  is  evaluated,  the  search  could  be  stopped  when  a  desired  quality  level  is  reached, 
rather  than  exhaustively  evaluating  the  experimental  design  points.  With  both  ele¬ 
ments,  the  search  could  be  conducted  as  parsimoniously  as  possible. 

6. 2. 3. 5  Automated  Decision  Agent.  Section  3. 2. 1.1  discusses  the  ex¬ 
perimental  design  built  to  investigate  a  range  of  values  of  aspiration  and  reservation 
levels.  If,  instead,  an  automated  decision  agent  (see  Sections  2. 1.1. 2  and  2.3.1)  could 
be  developed,  it  may  provide  better  insight  to  the  decision  maker.  In  fact,  even  the 
decision  strategies  of  a  decision  maker  could  be  investigated,  e.g.,  conservative  versus 
daring  decision  strategies. 

6.2.4  Software  Changes. 

6.2.4. 1  Code  to  Capture  Metrics.  Though  observations  have  been 
made  about  solution  quality  and  computational  efficiency,  because  existing  solvers 
were  used,  some  metrics  were  not  captured  and  could  not  be  reported.  The  following 
changes  to  the  software  are  suggested. 

1.  As  discussed  in  Section  4.4.1,  in  MGPS-RS,  the  difference  between  the  response 
standard  deviation  S  and  the  indifference  zone  parameter  is  used  in  the  ter¬ 
mination  criteria  for  the  algorithm;  however,  this  data  is  not  reported  at  the 
end  of  the  subproblem.  (Response  standard  deviation  in  this  case  is  that  of 
the  scalarized  objective  function,  not  the  individual  objective  functions.)  It  is 
suggested  that  the  response  standard  deviation  at  the  end  of  each  subproblem 
be  captured.  Though  not  a  direct  measure  of  the  variance  in  each  objective 
function,  it  does  provide  insight  into  the  variability  of  the  approximate  Pareto 
set. 

2.  Though  computational  efficiency  of  the  algorithm  was  discussed  in  Section  5.4, 
some  metrics  typically  associated  with  computational  efficiency  were  not  cap- 
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tured  by  the  existing  software  and  conld  not  be  reported.  Snch  measnres  wonld 
be  helpfnl  in  fnrther  assessing  the  algorithm’s  potential  nse  on  compntationally 
expensive  problems.  Specihcally,  it  is  snggested  that  the  software  be  modihed 
to  captnre  CPU  (central  processing  nnit)  time  nsed  by  each  snbproblem,  overall 
CPU  time,  nnmber  of  iterations  in  each  snbproblem,  and  nnmber  of  objective 
function  evalnations. 

6. 2. 4. 2  User-Friendly  Software.  This  research  provided  the  algorith¬ 
mic  methodology  for  a  two-stage  solntion  process  and  prototype  code  for  specihc 
implementations  of  the  algorithm.  This  code  was  integrated  with  the  batch  mode  of 
NOMADm,  bnt  was  not  integrated  with  the  graphical  nser  interface  (GUI).  In  fntnre 
research,  it  is  snggested  that  the  code  controlling  the  aspiration/reservation  level  anal¬ 
ysis  be  integrated  with  NOMADm  into  a  single  graphical  nser  interface.  Additionally, 
this  integrated  code  shonld  then  be  compiled  nsing  MATLAB®’s  VB.net®s  bnilder 
toolbox  so  that  it  can  be  execnted  by  nsers  withont  having  to  have  a  MATLAB® 
license. 
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Appendix  A.  Code — MOM  ADS  for  the  Dias  T2  Problem 


1.1  File  name:  DiasGamma2RuntheBatchFile.m 

SetUpProblem 

global  numobjectives ; 

global  aspire; 

global  reservation; 

global  utopia; 

global  nadir; 

global  objfunctionstruc ; 

global  xvals; 

“/oDefine  problem  files 

MyProb .problemPath  =  fullf ile (matlabroot , . . . 

’work’ , ’examples’ , ’theprob’ , ’ContBOlcky’ , ’DiasGamma2’ ) ; 

MyProb. F  =  ’DiasGamma2’ ;  %  functions  file 

MyProb. 0  =  ’DiasGamma2_0mega’ ;  °/o  linear  constraints  file 

MyProb. X  =  ’DiasGamma2_X’ ;  %  closed  constraints  file 

MyProb. 1  =  ’DiasGamma2_xO’ ;  °/o  initial  points  file 

MyProb. N  =  ’DiasGamma2_N’ ;  %  discrete  neighbor  file 

MyProb. P  =  ’DiasGamma2_Param’ ;  %  parameter  file 

MyProb. C  =  ’DiasGamma2_Cache .mat ’ ;  %  previously  created  Cache  file 

MyProb. S  =  ’DiasGamma2_Session.mat’ ;  %  previously  created  Session  file 

MyProb. H  =  ’DiasGamma2_History.txt’;  %  iteration  history  text  file 

MyProb. D  =  ’DiasGamma2_Debug.txt’;  %  debug  log  file 

MyProb. fType  =  ’M’ ;  %  type  of  functions  file  {M,F,C} 

MyProb. nc  =0;  %  number  of  nonlinear  constraints 

°/oOther  variables 

numobjectives  =  2;  utopia(l)  =  0;  utopia(2)  =  0;  nadir(l)  =  1; 
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nadir (2)  =  10; 

“/orange  checked 
arange= [.1  .8;  .05  1.]; 
rrange=[.81  1.2;  1 . 1  8]  ; 

numaspire=5;  nuiiireservation=5;  nuinreplications=5; 

“/Types  of  designs  ’fullfactorial’ ,  ’ centralcomposite ’  ,  ’boxbehnken’ 

“/You  need  the  statistics  toolbox  for  this  to  work 

“/Need  file  generatepoints .m  in  your  workpath 

TestPoints  =  generatepoints( ’fullfactorial’ ,numaspire, 

numreservation,  arange,  rrange,  numobjectives) ;  MyResult  =[]; 

objfunctionstruc  =  [] ;  bogus=0;  xvals  =  [] ;  for 

Mym=l : size (TestPoints , 1) 

for  Myn=l : numobjectives 

aspire(Myn)  =  TestPoints (Mym,Myn) ; 

reservation(Myn)  =  TestPoints (Mym,numobjectives+Myn) ; 

end 

if  (mod(Mym, 50)==0) 

Mym 

end 

for  1  =  1 rnumreplications 

temp  =  size (MyResult , 1)+1 ; 

“/This  is  where  the  MADS  or  GPS  solver  is  called. 

“/Need  the  MOmads.m  file  installed  in  the  workpath. 

“/For  MOmads,  you  need  the 
“/SetUpProblem.m  and  nomad. m  files  also. 

tempi  =  MOmads (MyProb, Options) ; 
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temp2=templ .par  am; 
if  temp  ==  1 
AddYes  =  1; 

else 

AddYes  =  ClieckPareto(temp2,  MyResult,  numobjectives) 

end 

if  AddYes  ==  1 

MyResult (temp, 1) .X  =  tempi. x; 

MyResult (temp, 1) . f  =  tempi. f; 

MyResult (temp, 1) . c  =  tempi. c; 

MyResult (temp , 1) . aspire=aspire ; 

MyResult (temp , 1) . reservation=reservation; 

MyResult (temp, 1) . 0bjective=temp2; 

else 

bogus  =  bogus+1; 

end 

end 

end 

°/oThis  section  was  for  graphing  the  output.  Uncomment  as  needed. 

°/o  for  m=l :  size  (MyResult,  1) 

°/o  MyX(m)=MyResult(m,  1)  .x(l) ; 

°/o  MyY(m)=MyResult(m,  1)  .x(2) ; 

°/o  end 

°/o  scatter (MyX,MyY, 3, ’filled’ ) 

°/o  clear  MyX  MyY ; 

°/„  figure 
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%  hold 

°/o  for  m=l :  size(MyResult,  1) 

°/o  MyX(m)=MyResult (m,  1)  . Objective (1) 

°/o  MyY(m)=MyResult (m,  1)  . Objective (2) 

°/o  end 

°/o  scatter (MyXjMyY, 3, ’filled’ ) 
y.  clear  MyX  MyY; 

°/„  figure 
y.  hold 

°/o  clear  MyX  MyY ; 

°/„  figure 
y.  hold 

°/o  for  m=l :  size(MyResult,  1) 

°/o  MyX (m)=MyRe suit (m,  1)  . Objective (2) 

°/o  MyY(m)=MyResult (m,  1)  . Objective (3) 

°/o  scatter(MyX,MyY,3, ’filled’ ) 

°/o  end 

y.  clear  MyX  MyY; 
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1.2  File  name:  DiasGamma2.m 


function  [fx,cx]  =  DiasGamma2(x) ; 
global  numobjectives ; 
global  aspire; 
global  reservation; 
global  utopia; 
global  nadir; 
global  objfunctionstruc ; 
global  xvals; 

°/of(l)  and  f(2)  are  the  objective  functions, 
f (l)=x(l)+.01*rand;  gl=sum(x)-x(l) ;  g2=l+(9/29) *gl ; 
f (2)=g2*(l-(f (l)/g2)~2)+0.01*rand; 

“/oThis  saves  the  original  function  values  to  the  variable  param. 

“/oUsed  in  the  batch  file  for  the  filter. 

Param . param=f ; 

°/oSet  epsilon  and  alpha  and  beta  and  w 
epsilon  =  5; 

“/oThis  code  puts  the  objectives  into  the  achievement 
°/o  scalarization  function, 
for  nl  =  1 : numobjectives 

w(nl)=  1/ (reservation(nl) -aspire (nl)) ; 
if  (aspire (nl)  ~=  utopia(nl)) 

alpha (nl) =0 . 1* (reservat ion (nl) -aspire (nl) ) / (aspire (nl) -utopia (nl) ) ; 
else 

alpha(nl)=0. l*(reservation(nl) -aspire (nl))/ (0.0000001) ; 
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end  if  (aspire (nl)  ~=  nadir (nl)) 

beta(nl)=-10* (reservation(nl) -aspire (nl) )/ (aspire (nl)-nadir (nl) ) ; 
else 

beta(nl)=-10* (reservation(nl) -aspire (nl) )/ (0 . 0000001) ; 
end  end  for  n2  =  1 :numobjectives 
if  f(n2)  <  aspire (n2) 

u(n2)=alpha(n2)*w(n2)* (aspire (n2)-f (n2) )+l ; 
elseif  f (n2)<=reservation(n2) 

u(n2)=w(n2) * (aspire (n2)-f (n2) )+l ; 

else 

u(n2)=beta(n2)*w(n2)*(reservation(n2)-f (n2)) ; 

end 

end  temp  =  [u(l)  u(2)] ; 

°/oThis  is  the  combined  objective  function, 
fx  =  -(min(temp)+epsilon*sum(temp) ) ; 

°/oThis  would  be  used  if  you  had  non-linear  constraints. 
cx=  []  ; 

°/oDon^t  delete  this  line  or  the  param  variable  above  won’t  save. 
setappdata(0, ’PARAM’ , Param) ; 

return 
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1.3  File  name:  DiasGamma2_Omega.m 

°/oCanoeDW_Omega:  User-supplied  funct.  for  defining  Omega,  based  on  p. 

function  [A,l,u,plist]  =  DiasGamma2_0mega(n) ;  A  =  eye(n);  1  =  [0;  0; 

0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  ... 

0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0]; 
u  =  n-  1-  1-  1-  1-  1-  1-  1-  1-  1-  1-  1-  1-  1-  1-  1-  1- 
1-  1-  1-  1-  1-  1-  1-  1-  1-  1-  1-  1-  ll- 

X,  X,  X,  X,  X,  X,  X,  X,  X,  X,  X,  X,  XJ, 

return 

1.4  File  name:  DiasGamma2-xO 

function  iterate  =  DiasGamma2_xO; 

°/oThis  sets  the  initial  points. 

iterate (l).x  =  [0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  ... 

0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0;  0]; 

iterate (l).p  =  {};  return; 
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1.5  File  name:  SetupProblem.m 

y.  DO  NOT  MODIFY  THIS  SECTION 

°/o  Set  Options  to  their  default  values 

cic  clear  variables  Defaults  =  mads_defaults( ’Truth’ ) ;  Options 
Defaults . Options ; 

“/.Options  MAKE  CHANGES  HERE  FOR  SEARCH  TYPE,  ETC. 

“/.  Specify  Choices  for  SEARCH 
Options .nSearches  =  2; 

Options . Search(l)  .type  =  ’LHS’;"/.  For  choices,  see  mads_def aults 

Options . Search(l)  .niter  =  1;  “/.  Number  of  iter  for  Search  #1 

Options . Search(l)  .nPoints  =  8;  “/.  Number  of  poll  or  sample  points 

Options .  Search(l)  .  sfile  =  ’’;  “/.  filename  must  include  full  path 

Options . Search(l)  . file  =  “/.  filename  must  include  full  path 

Options . Search(l)  . local  =  0;  “/.  flag  to  turn  on  trust  region 

Options . Search(l)  .merit  =  0;  “/.  flag  to  penalize  clustered  data 

Options . Search (1)  .param  =  0;  “/.  flag  to  constr.  surr.  on  param 

Options . Search(l) . recal  =  ... 

strncmp (Options . Search(l) . type (max (1 , end- 1) : end) ,  ... 

’NW’ ,2)  II  strncmp (Options . Search(l) . type (max (1 , end-3) ;end) , ’DACE’ ,4) ; 
Options . Search(2)  .type  =  ’None’;  “/.  For  choices,  see  mads_def aults 

Options . Search(2)  .niter  =  10;  “/.  Number  of  iter  for  Search  #2 

Options . Search(2)  .nPoints  =  1;  “/.  Number  of  poll  or  sample  points 

Options .  Search(2)  .  sfile  =  ’regpoIyO’ ;“/.  filename  must  include  full  path 
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Options . Search(2)  . file  =  °/o  filename  must  include  full  path 

Options . Search(2)  . local  =  0;  °L  flag  to  turn  on  trust  region 

Options . Search(2) .merit  =  0;  %  flag  to  penalize  clustered  data 

Options . Search(2) .recal  =  ... 

strncmp(Options.Search(l) .type(max(l, end-1) :end) ,  ’NW’,2)  ||  ... 

strncmp (Opt ions . Search (1) . type (max (1, end-3) :end) ,  ^DACE’ ,4) ; 

Options . SurOptimizer  =  ’mads’;  Options .mvplSurrogate  =  1; 

0ptions.dace(2) .reg  =  ’regpolyO’;  Options . dace (2) . corr 

’correxp’;  Options .dace (2) .theta  =  1;  Options. dace (2) .lower 
=  0.01;  Options. dace (2) .upper  =  1000;  Options .dace (2) . isotropic 
=  0; 

°/o  Specify  Choices  for  POLL 


Options . pollStrategy 

= 

’MADS_2n’;  “/. 

For  choices,  see  : 

mads_def aults 

Options . pollOrder 

= 

’Consecutive 

’;  7  For  choices, 

see  defaults 

Options . pollCenter 

= 

0;  7o 

Poll  around  n-th 

filter  point 

Options . pollComplete 

= 

0;  7 

Flag  for  complete  polling 

Options . NPollComplete 

= 

0;  7 

Flag  for  complete 

neigh,  polling 

Options . EPollComplete 

= 

0;  7 

Flag  for  complete  extend  polling 

°/o  Specify  Termination 

Criteria 

Options . Term . delta 

= 

le-4;  7 

minimum  mesh  size 

Options .Term. niter 

= 

Inf ;  7 

maximum  number  of 

iterations 

Options .Term. nFunc 

= 

50000 ;  7 

maximum  number  of 

function  evals 

Options .Term. time 

= 

Inf ;  7 

maximum  CPU  time 

Options .Term. nFails 

= 

Inf ;  7 

max  number  of  consec.  Poll  fails 

°/o  Choices  for  Mesh  Control 
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Options . deltaO  =  1; 

Options . de It aMax  =  1; 

Options .meshRefine  =  0.5; 

Options .meshCoarsen  =  2.0; 


°/o  Choices  for  Filter  management  (for 
Options.hmin  =  le-8; 

Opt ions. hmax  =  1.0; 

%  Choices  for  EXTENDED  POLL  (for  MVP 
Options . ePollTriggerF  =  0.01; 

Options . ePollTriggerH  =  0.01; 

°/o  MADS  flag  parameter  values 


Options . loadCache  =  0; 

Options . countCache  =  0; 

Options . runStochastic  =  1; 

Options . scale  =  2; 

Options .useFilter  =  1; 

Options . degeneracyScheme  =  ’random’; 
Options . removeRedundancy  =  1; 

Options . computeGrad  =  0; 

Options . saveHistory  =  0; 

Options .plotHistory  =  0; 

Options .plotFilter  =  0; 

Options .plotColor  =  ’k’; 

Options . debug  =  0; 

Options . showFilterPlot  =  0; 


°/o  initial  mesh  size 

°/o  bound  on  how  coarse  mesh  can  get 

°/o  mesh  refinement  factor 

°/o  mesh  coarsening  factor 

problems  with  nl  constraints) 

°/o  minimum  infeasible  point  h-value 
°/o  maximum  h-value  of  a  filter  point 

problems) 

°/o  f-value  Extended  Poll  trigger 
°/o  h-value  Extended  Poll  trigger 

°/o  load  pre-existing  Cache  file 
°/o  count  Cache  points  as  func.  calls 
%  runs  problem  as  a  stoch.  problem 
°/o  scale  directions  using  log  base 
%  filter (0=none , l=multi-pt ,2=2-pt) 

°/o  scheme  for  degenerate  constraints 
°/o  discard  redundant  lin.  constr. 

°/o  compute  gradient ,  if  available 

°/o  saves  MADS  perform,  to  text  file 

°/o  plot  MADS  performance 

°/o  plot  the  filter  real-time 

°/o  color  of  history  plot 

°/o  turn  on  status  mess  for  debugging 

°/o  turn  off  the  filter  plot 
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Opt ions. Sur.Term.de It a  =  0.01; 


°/o  surrogate  minimum  mesh  size 


°/o  Set  up  figure  handles  for  real-time  plots 


if  (Options .plotHistory  ==  2) 
f igure ; 

Options .hplothandle  =  gca; 

end 
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1.6  File  name:  generatepoints.m 

function  [TestPoints]  =  generatepoints(designof experiments, 
numasp , . . . 

numres,  arange,  rrange,  numobj) 
switch  lower (designof experiments) 
case  ’fullfactorial^ 
gen  =  []  ; 
for  ar=l:2 

for  n=l: numobj 
if  ar==l 

gen  =  [gen  numasp]  ; 

else 

gen  =  [gen  numres] ; 

end 

end 

end 

myf act=fullf act (gen) ; 

testpt=myf act ; 

for  m=l : size (myf act , 1) 

for  1=1 : size (myf act , 2) 
if  l<=numobj 

testpt (m,l)=arange(l , l)+(myf act (m, 1)-1) * . . . 

( (arange (1 , 2) -arange (1 , 1) )/ (numasp-1) ) ; 

else 

testpt(m,l)=rrange(l-numobj , l)+(myfact(m,l)-l) . . . 
* ( (rrange (1 -numobj , 2) - . . . 
rrange (1 -numobj , 1) )/ (numres-1) ) ; 

end 

end 
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end 


TestPoints=testpt ; 
case  ’ central composite’ 

myf act=ccdesign(2*numobj ) ; 

testpt=myf act ; 

for  m=l : size (myf act , 1) 

for  1=1 : size (myf act , 2) 
if  l<=numobj 

testpt (m,l)=0 . 5* ( (aranged ,2)+arange (1,1))  +  -  • 
(myf  act  (m,  1)  *  (aranged  ,2) “grange  (1,1)))); 

else 

testpt (m,l)=0. 5* ((rranged-numobj  ,2)  +  .  .  . 
rranged-numobj  ,  1) )  + (myf  act  (m,  1)  * .  .  . 
(rrange (1-numobj , 2) -rrange (1-numobj , 1) ) ) ) 

end 

end 

end 

TestPoints=testpt ; 

case  ’boxbehnken’ 

myf act=bbdesign(2*numobj ) ; 

testpt=myf act ; 

for  m=l : size (myf act , 1) 

for  1=1 : size (myf act , 2) 
if  l<=numobj 

testpt (m,l)=0 . 5* ( (aranged ,2) +arange (1,1))  +  .  . 
(myf  act  (m,  1)  *  (aranged  ,2) -arange  (1,1)))); 

else 

testpt (m,l)=0. 5* ((rrainge (1-numobj  ,2)  +  .  .  . 
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rrange(l-numobj , 1) )+(myf act (m, 1) * . . . 
(rrange (1-numobj , 2) -rrange (1-numobj , 1) ) ) ) 

end 

end 

end 

TestPoints=testpt ; 
case  ’singlepoint’ 
testpt=  []  ; 
for  n=l;2*numobj 
if  n<=numobj 

testpt= [testpt  arange(n, 1)] ; 

else 

testpt= [testpt  rrangeCn, 1)]  ; 

end 

end 

TestPoints=testpt ; 
otherwise 

dispC’I  need  a  valid  type  of  design  of  experiments.’) 

end 
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1.7  File  name:  CheekPareto.m 

function  YesNo  =  CheckPareto(functionvalues,  MyResult, 
numobj actives)  thisone=[] ; 

for  m=l :numobj active s 

thisone  =  [thisone  functionvaluesd ,m)]  ; 
end  for  n=l : size (MyResult , 1) 
that  one  =  []  ; 
for  m=l :numobj actives 

thatone  =  [thatone  MyResult (n, 1) -Obj active (m)] ; 

end 

MyCheck=sum(su]n(thatone  <  thisone)); 
if  MyCheck  ==  numobj actives 
YesNo=0 ; 
break 

end 

YesNo=l; 

end 

1.8  File  name:  MOmads.m 

function  [BestF,BestI ,RunStats ,RunSet]  =  MOmads(MyProb, Options) 
°/oMADS_BATCH  Sets  up  and  runs  the  MADS  algorithm  without  a  GUI . 

y. 

problemPath  =  MyProb.problemPath; 

Problem. File .F  =  MyProb.F;  %  functions  file 

Problem. File . 0  =  MyProb.O;  °/o  linear  constraints  file 

Problem. File .X  =  MyProb.X;  °/o  closed  constraints  file 

Problem. File . I  =  MyProb.I;  %  initial  points  file 
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Problem. File .N  =  MyProb.N; 
Problem. File .P  =  MyProb.P; 
Problem. File. C  =  MyProb.C; 
Problem. File. S  =  MyProb.S; 
Problem. File .H  =  MyProb.H; 
Problem. File .D  =  MyProb.D; 
Problem . f Type  =  MyProb . f Type ; 
Problem. no  =  MyProb. no; 
Problem. nameCache  =  ’CACHE’; 


°/o  discrete  neighbor  file  (MVP  only) 

°/o  parameter  file 

°/o  previously  created  Cache  file 

°/o  previously  created  Session  file 

°/„  iteration  history  text  file 

°/o  debug  log  file 

°/o  type  of  functions  file  {M,F,C} 

°/o  number  of  nonlinear  constraints 
Problem. typeProblem  =  ’TRUTH’; 


y.  DO  NOT  MODIFY  AFTER  THIS  POINT 

if  (Problem. nc  ==  0) 

Options .plotFilter  =  0; 

end  if  (Options .plotFilter)  &&  (Options . showFilterPlot  ==  1) 
f igure ; 

Options . fplothandle  =  gca; 

end 

%  Set  the  path,  and  load  any  user-provided  problem  parameters 
cwd  =  pwd;  “/ocreate  variable  to  remember  current  directory 

cd(problemPath) ;  Zchange  focus  to  problem  directory 
if  (exist (Problem. File . P, ’file ’ )  ==  2) 

Problem. Param  =  feval (Problem. File .P) ; 
setappdata(0, ’PARAM’ , Problem. Param) ; 

end 
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°/o  Get  the  initial  iterates  and  call  the  optimizer 
if  isf  ield(Problem, ’ParamO  &&  isf  ield(Problem.Param,  ’ iterateOO 
iterateO  =  Problem. Param. iterateO; 
else 

iterateO  =  feval(Problem.File . I) ; 

end 

[BestF,Bestl,RunStats,RunSet]  =  mads (Problem, iterateO, Options) ; 


cd(cwd) ;  “/ochange  focus  back  to  original  directory 

return 
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